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Abstract 



Using the advective Cahn-Hilliard equation as a model, we illuminate the role of advection in 
phase-separating binary liquids. The advecting velocity is either prescribed, or is determined by 
an evolution equation that accounts for the feedback of concentration gradients into the flow. 

After obtaining some general results about the existence and regularity of solutions to the model 
equation, we focus on two specific cases: advection by a chaotic flow, and coupled Navier-Stokes 
Cahn-Hilliard equations in a thin geometry. 

By numerically simulating chaotic flow, we show that it is possible to overwhelm the segregation 
by vigorous stirring, and to create a homogeneous state. We analyze the mixing properties of the 
model: by measuring fluctuations of the concentration away from its mean value, we find a priori 
bounds on the amount of homogenization achievable. 

We discuss the Navier-Stokes Cahn-HiUiard equations and derive a thin-film version of these 
equations. We examine the dynamical coupling of the concentration and velocity (backreaction). 
To study long-time behaviour, we regularize the equations with a Van der Waals potential. We 
obtain existence and regularity results for the thin-film equations; the analysis also provides a 
nonzero lower bound for the film's height, which prevents rupture. 

We carry out numerical simulations of the thin-film equations and, by comparing the results with 
experiments on polymer blends, we show that our model captures the qualitative features of real 
binary liquids. The outcome of the phase separation depends strongly on the backreaction, which 

we demonstrate by applying a shear stress at the film's surface. When the backreaction is small, the 
domain boundaries align with the direction of the stress, while for larger backreaction strengths, 
the domains align in the perpendicular direction. 

Lastly, we compare and contrast the Cahn-Hilliard equation with other models of aggregation; this 
leads us to investigate the orientational Holm-Putkaradze model. We demonstrate the emergence 
of singular solutions in this system, which we interpret as the formation of magnetic particles. 
Using elementary dynamical systems arguments, we classify the interactions of these particles. 
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1.1 Spinodal decomposition with wetting. Subfigure (a) provides a schematic represen- 
tation of the spinodal wave at the lower horizontal wall. This structure forms when 
one of the binary fluid components wets the wall. Layers of alternating demixed 
regions propagate into the bulk until the wall effects are negligible. This effect is 
suppressed in thin films. Subfigure (b) shows the completely wet thin film, in which 
a single binary fluid component wets the horizontal bounding walls. The film is too 
thin to support a spinodal wave. Subfigure (c) shows the partially wet thin film, in 
which both binary fluid components wet the horizontal bounding walls, and there is 
no vertical domain structure. In this report we shall focus on the partially wet case. [19] 

3.1 (a) Lyapunov exponent 70 and its small- and large-amplitude forms; (b) Comparison 
between the Lyapunov exponent and 72; (c) Comparison between the Lyapunov 
exponent for the sine map with and without the randomization of phases [53] 



3.2 (a) Concentration, at t = 30,000; (b) Structure function kfs {k/ki,t) in the scaling 
regime t > 20, 000; (c) Time dependence of ki, with power law exponent close to the 
LS exponent 1/3; (d) Free energy, F, exhibiting the time dependence F ~ 
close to the LS exponent 



3.3 (a) F vs t for (from bottom) a = 0,0.1,0.3,0.5,0.7,1.0; (b) ki vs t for the same 
values of a; (c) Tracer variance for (from bottom) a = 1.0,0.7,0.5,0.3,0.1,0. . . [57] 

3.4 The steady-state concentration field after t = 30, 000 timesteps, for various stirring 
amplitudes a [58] 

3.5 (a) The plot of log(cr^/F) against log A suggests the scaling law (a'^/F) ~ X^^^^ for 
small stirring amplitudes while for large amplitudes this proxy for the bubble radius 
decays exponentially with increasing Lyapunov exponent, (b) Graph of log(/ci) 
against log A, with approximate power law behaviour (ki) ~ A°'^^. The exponent 
is not the same over the entire range, however, (c) The graph of log(-F) exhibits 
a much cleaner power law behaviour {F) ~ X^'"^^. (d) The Batchelor wavenumber 

fcfe = [|Vc|Vc2] also possesses a clear power law behaviour {kf,) ~ A°'^^ [59] 
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3.6 Normalized PDF of concentration in the steady state, (a), (b) Segregation-dominated 
flow (a = 0.1,0.3); (c), (d) Crossover to quasi-diffusive regime (a = 0.5,0.7); (e) 
Quasi- diffusive regime (a = 1.0); (f) Semilog plot of PDF for a = 1.0. Note Gaus- 
sian core and exponential decay down to the cutoff at |c| < 1 
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the homoclinic case, the zeros of lie at the turning points Ci and C2 EH 
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the attractor (co,0); (b) Profile of the solution and its derivative; (c) The full two- 
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3.11 (a) Front profile for Xj/D = 0.001,0.01,0.1, 1, 10, 100, 1000. The front steepens as 
X'-f/D increases; (b) Front width as a function of the parameter X'y/D: for large 
parameter values, the front width scales hyperdiffusively as (X'j/D)^^^ 



3.12 (a) The plot of log((7^/F) against log A suggests the scaling law (a'^/F) ~ X~^^'^ 
for small stirring amplitudes, while for large amplitudes this proxy for the bubble 
radius decays exponentially with increasing Lyapunov exponent, (b) The graph of 
log(-F) exhibits power law behaviour (F) ~ A°'^^. The predominance of kinetic free 
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the dependence of m™™ on the velocity amplitude Vq. In (a), the scale of the source 
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5.1 Definition sketch for the long- wavelength approximation in two dimensions. The 
approximation relies on the smallness of the parameter 6 = h^/X 



5.2 A sketch of the regularized function (s) obtained in Eq. (5.21) 



6.1 Equilibrium solutions of the thin-film equations in one dimension, obtained from 
solving a boundary- value problem with parameter values C = C'^\A\ = 1 and r = 
0.1,1,10,50. In (a) the valley deepens with increasing r although the film never 
ruptures, while in (b) the front steepens with increasing r. In (c) we plot the forces 
Fcap and Fydw for C = = |v4| = 1 and r = 50 to verify that they have opposite 

sign 11161 
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and (b) show the dip magnitude's dependence on the parameter C^|A|; Subfigures 
(c) and (d) give the behaviour of the dip magnitude as a function of r. In both 
cases, the approximate relation h~^^^ ~ r/Cn|y4| is discernable 11171 

6.3 A typical plot of 5 C) for Fq = Fi = |^and C = 1. This theoretical lower bound 



has a different shape from those in Fig. 6.2, which suggests that while 



plays an important role in the analysis of the model equations, it does not capture 

the physics of film thinning 11181 

6.4 (a) A numerical solution of the SCH equations in one dimension, with C = —A = 1 
and r = 1 (large backreaction); (b) The evolution of the free-energy functional. The 
free energy decays monotonically, in agreement with the theory, except at late time, 
where there is a slight loss of monotonicity due to numerical error 11191 

6.5 Log-log plot indicating the time dependence of the typical domain size R{t). This 
dependence is close to the Lifshitz-Slyozov growth law, R ~ t^^^ 11211 

6.6 (a) Concentration field at t = 5,000, r = 0.01, C = 1, and A = -10; (b) The 
height field for the same parameter values. Note that the height field is slaved to 

the concentration field 11221 
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6.7 Domain structure at t = 1000, t = 3750, t = 7500, and t = 15000 for C = 1, 
A = —10, and r = 0. The surface tension gradient is parallel to the arrow and 

r = Fq sin (kx), Tq = 20 and k = 4fco- For r = 0, the domains align along the arrow. 11231 

6.8 Domain structure at t = 1000, t = 3750, t = 7500, and t = 15000 for C = 1, 
A = —10, and r = |. The surface tension gradient is parallel to the arrow and 
r = Fq sin (kx), Tq = 20 and k = 4fco- The system reaches a steady state where the 
domains align perpendicular to the arrow 11231 

6.9 The free-surface height for r = and t = 15000 aligns with the applied surface 
tension. The height field at t = 15000 for r = | is similar 11241 

6.10 Time dependence of the dominant spatial scales kx and ky for (a) r = 0, where the 
secular behaviour of k^ and ky has a small drift whose rate we do not report; (b) 

r = |, where k^ ^ 2 and ky ^ 11251 



7.1 The initial data for the magnetization equation (7.3). This initialization is obtained 



by allowing the orientation angles of the magnetization vector to vary sinusoidally 



in space, as in Eq. (7.10). Here the wavenumber of the variation is equal to the 



fundamental wavenumber 2tt/L 11311 

7.2 Numerical simulations of Case (1), the Landau-Lifshitz-Gilbert equation in the 
over-damped limit. In this case, the magnetization decays to a constant state. 
Subfigures (a) and (b) show the magnetization at times t = 0.03 and t = 0.15 
respectively; (c) is the energy functional, which exhibits exponential decay after 
some transience. The final orientation is (0, 6*) = (tt, 7r/2) 11321 

7.3 Numerical simulations of Case (2), the nonlocal Gilbert equation with with a < f3. 
In this case, the energy increases to a constant value, and the magnetization becomes 
constant. Subfigures (a) and (b) show the magnetization at times t = 8 and t = 40; 

(c) is the energy functional. The final orientation is (0, 6') = (vr, 7r/2) 11321 

7.4 Numerical simulations of Case (3), the nonlocal Gilbert equation with a > (3. In 
this case, the energy decreases indefinitely, and the magnetization vector develops 
finer and finer scales. Subfigures (a), (b), and (c) show the magnetization at time 
t = 10000; (d) is the energy functional, which decreases in time as a power law at 
late times. Subfigure (e) shows the power spectrum of m^; the integer index n labels 
the spatial scales: if /c„ is a wavenumber, then the corresponding integer label is 

n = knL/2TT [133] 

7.5 The first route to instability. Subfigure (a) shows the growth rate cxi for < Pm < 
a^, with negativity indicating a stable equilibrium; (b) gives the growth rate ai for 
oiip < am < Pm, with positivity indicating an unstable equilibrium. We have set 

\mo\ =i!o = l [136] 
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7.6 The second route to instability. Subfigure (a) shows the growth rate (J2 for am < Pm, 
with negativity indicating a stable equilibrium; (b) gives the growth rate a2 for am > 

j3m, with positivity indicating an unstable equilibrium. We have set |mo| = ipo = 11371 

7.7 The nuUcline dx/dt = of the two-clumpon dynamical system with am < 01 p. 
The region contained inside the dotted lines ?/ = ±1 gives the allowed values of 
the dynamical variables {x,y). Subfigure (a) shows the case when Pm < C(m- The 
stable equilibria of the system are (x, y) = (±(i, 1) and the line a; = 0. All initial 
conditions flow into one of these equilibrium states; subfigure (b) shows the case 
when am < Pm- Initial conditions confined to the line y = 1 flow into the fixed 
point (±(i, 1), while all other initial conditions fiow into the line x = 11391 

7.8 The nullcline dx/dt = of the two-clumpon dynamical system with ap < am- The 
region contained inside the dotted lines y = ±1 gives the allowed values of the 
dynamical variables {x,y). Subfigure (a) shows the case when Pm < Cim- The lines 
X = and x = ±00 form the stable equilibria of the system. All initial conditions 
flow into one of these states; subflgure (b) shows the case when am < Pm- Initial 
conditions conflned to the line y = 1 flow into the flxed points (0, 1) and (±00, 1), 
while all other initial conditions flow into the line a; = 11401 

7.9 Evolution of sinusoidally- varying initial conditions for the density and magnetiza- 



tion, as in Eq. (7.22). Subflgure (a) shows the evolution of * ip for t G [0, 0.15], 
by which time the initial data have formed two clumpons; (b) shows the evolution of 
fj^x- The proflles of fiy and /i^ are similar. Note that the peaks in the density proflle 
correspond to the troughs in the magnetization proflle. This agrees with the linear 
stability analysis, wherein disturbances in the density give rise to disturbances in 
the magnetization 11421 

7.10 Evolution of sinusoidally- varying initial conditions for the density and magnetiza- 



tion, as in Eq. (7.22). Subflgure (a) shows the system velocity V given in Eq. (7.12), 
just before the singularity time; (b) shows the magnetization Hm at the same time. 
The density maxima emerge at the locations where the convergence of —V (flow 
into X = and x = ±L/2) occurs, and the magnetization develops extrema there. . 11431 

7.11 Evolution of a flat magnetization fleld and a sinusoidally- varying density, as in 



Eq. (7.23). Subflgure (a) shows the evolution of H^*ip for t G [0.5, 1]; (b) shows the 
evolution of ^x- The proflles of /ly and fiz are similar. At t = 0.5, the initial data 
have formed eight equally spaced, identical clumpons, corresponding to the eight 
density maxima in the initial conflguration. By impulsively shifting the clumpon 
at X = by a small amount, the equilibrium is disrupted and the clumpons merge 
repeatedly until only one clumpon remains 11441 

7.12 Evolution of a flat magnetization fleld and a sinusoidally- varying density, as in 



Eq. (7.23). Subflgure (a) shows the system velocity V given in Eq. (7.12) just 
before the singularity time; (b) gives the magnetization Hm at the same time. The 
density maxima emerge at the locations where the convergence of —V occurs, and 
the magnetization develops extrema there 11451 
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7.13 Evolution of a flat magnetization field and a sinusoidally- varying density, as in 
Eq. (7.23). Shown is the velocity profile for t G [0.5, 1]; the system velocity is given 
by Eq. (7.12). The velocity —V flows into each density maximum, concentrating 
matter at isolated points and precipitating the formation of eight equally-spaced 
identical clumpons. On a periodic domain, such an arrangement is an equilibrium 
state, although it is unstable. Thus, by impulsively shifting the clumpon at a; = 
by a small amount, we force the clumpons to collapse into larger clumpons, until 
only a single clumpon remains 
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Chapter 1 
Introduction 



Overview, review, preview 

In her column in The Guardian, Maureen Lipman has depicted the routine of making home-made 
mayonnaise: 

One day she... gave me a demonstration on how to make mayonnaise. I had no 
idea it was so technical... She whisked the mustard with one yolk for a few minutes, 
then started dribbling in the oil. As soon as any separation appeared she whisked even 
faster and continued whisking and oiling for long enough to make my wrist hurt, let 
alone hers. It was riveting, like watching an old master mixing his ochres with his 
burnt siennas. pL| 

In this report we shall study this process in a more abstract setting, together with the parallel 
problem of phase separation in thin layers. We wish to characterize the competing effects of ho- 
mogenization and phase separation in a stirred binary fluid. The advective Cahn-Hilliard equation 
describes precisely these mechanisms and this will be the subject of our report. In this chapter 
we shall review the literature pertaining to both stirred and unstirred phase separation under 
Cahn-Hilliard (CH) dynamics and then outline the study we have made. 

In a seminal paper, Cahn and Hilliard [2] introduced their eponymous equation to model the 
dynamics of phase separation. They pictured a binary alloy in a mixed state, cooling below a 
critical temperature. This state is unstable to small perturbations so that fluctuations cause the 
alloy to separate into bubbles (or domains) rich in one material or the other. A small transition 
layer separates the bubbles. Because the evolution of the concentration field is an order-parameter 
equation, this description is completely general. Thus, by coupling the phase-separation model to 
fluid equations, one has a broad description of binary fluid flow encompassing polymers, immiscible 
binary fluids with interfacial tension, glasses, and mayonnaise. 

In industrial applications, the components of a binary mixture often need to be mixed in a ho- 
mogeneous state [3j. In that case, the coarsening tendency of multiphase fluids is undesirable. In 
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this report we shall examine how a stirring flow affects this mixing process. We are interested in 
stirring not just to limit bubble growth, as is often effected with shear flows, but to break up the 
bubbles into a homogeneous mixture. We shall also examine the interplay between the concentra- 
tion gradients and the fluid velocity, and shall find that this coupling effect can play an important 
role in the control of phase separation 

We outline previous work for three increasing levels of complexity: CH in the absence of flow, 
passive CH with flow, and active CH with flow. 

Cahn-Hilliard fluids without flow: There is a substantial literature that treats of the CH equation 
without flow. For a review, see The main physical feature of this simpler model is the 
existence of a single lengthscale (the bubble or domain size) such that quantities of interest (e.g., 
the correlation function) are self-similar under this scale. This length grows in time as t^^^, a result 
seen in many numerical simulations (see, for example, Zhu et al. [S]). Of mathematical interest 
are the existence of a Lyapunov functional, and hence of a bounded solution [6], and the failure 
of the maximum principle owing to fourth-order derivatives in the equation for the concentration 

field [giH]. 

Passive Cahn-Hilliard fluids: Given an externally-imposed flow (e.g., a shear flow, turbulence, 
or a stirring motion), we have a passive tracer equation for a CH fluid. There have been several 
studies of the (active and passive) shear-driven CH fluid PHH]. In these studies, it is claimed that 
domain elongation persists in a direction that asymptotically aligns with the flow direction, while 
domain formation is arrested in a direction perpendicular to this axis. This seems to be confirmed 
in the experiments of Hashimoto pT] , although these results could be due to finite-size effects [12] , 
in which the arrest of domain growth is due to the limitation of the boundary on domain size. Of 
interest is the existence or otherwise of an arrest scale and the dependence of this scale on the 
system parameters. 

A more effective form of stirring involves chaotic flows, in which neighbouring particle trajectories 
separate exponentially in time, leading to chaotic advection [13]. The average rate of separation 
is the Lyapunov exponent, positive for such a flow [TH [T5] . One may also regard this exponent 
as the average rate of strain of the flow. By imposing a chaotic flow on the CH fluid, finite-size 
effects are easily overcome. This is because the shear direction changes in time, so that domain 
elongation in one particular direction, leading to an eventual steady state in which domains touch 
many boundaries, is no longer possible. Berthier et al. [16] investigated the passive stirring of the 
CH fluid and observed a coarsening arrest arising from a dynamical balance of surface tension and 
advection. The stirring by the chaotic flow destroys large-scale structures, and this is balanced 
by the domain-forming tendency of the CH equation. A balancing scale is identified with the 
equilibrium domain size. 

The importance of exponential stretching is amplified by the results of Lacasta et al. [T7], in 
which a passive turbulent flow endowed with a tunable amplitude, correlation length, and cor- 
relation time, causes both arrest and indefinite growth, in different limits. We explain this 
heuristically as follows: In one limit of turbulent flow (small Prandtl number), velocity corre- 
lations are short in range compared to bubble size, fluid particles diffuse [18], and bubble radii 
evolve as dR^/dt = 2k,cs + [Cahn-Hilliard contribution], where Rb is the bubble radius and Hes is 
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the effective tracer diffusivity. Because Kes > 0, the presence of diffusion fails to arrest bub- 
ble growth. In the limit of large Prandtl number however, the velocity field has long-range 
correlations and exhibits exponential separation of trajectories, so that the radius equation is 
dRf/dt = —2XRl + [Cahn-Hilliard contribution], where A is the average rate of strain, or Lya- 
punov exponent, of the flow. Thus, the exponential stretching due to stirring balances the bubble 
growth. 

Active Cahn-Hilliard fluids: At the highest level of complexity, one considers the reaction of the 
mixture on the flow by coupling the CH equation to the Navier-Stokes equations. This choice of 
coupling is not unique, although the different models agree in the case that interests us, namely 
density-matched fluids [121 [20] • In experiments involving turbulent binary fluids near the critical 
temperature, it is found that a turbulent flow suppresses phase separation and homogenizes the 
fluid [211 [22] . However, by cooling the fluid further and maintaining the same level of forcing, 
phase separation is achieved. This result is limited to near-critical fluids, although it does suggest 
the possibility of homogenization by stirring in other contexts. 

The most recent work on the CH fluid (Berti et al., [20j) focusses on the active CH tracer. They 
couple the CH concentration field to an externally-forced velocity field and observe coarsening 
arrest independently of finite-size effects. They identify a balance of local shears and surface 
tension as the mechanism for coarsening arrest. Thus, the chaotic advection discussed above is 
a sufficient condition for coarsening arrest, but not a necessary one since shear will suffice. Of 
interest again is the dependence of the arrest scale on the mean rate of shear. 

With this literature in mind, we study the case of a chaotic flow coupled passively to the CH 
equation. Our ultimate focus will be on the breakup of bubbles by stirring, and not just coarsening 
arrest. We believe this aspect of CH flow deserves a further and more complete study, as it is a 
regime of particular relevance to industry, for example, in droplet breakup or in the mixing of 
polymers [23] . 

Berthier et al. [16] used an alternating sine flow to study coarsening arrest. Because of its simplicity, 
the sine flow is a popular testbed for studying chaotic mixing [211 [SSI [261 [23 [28]. In contrast to 
Berthier et al. [TU], we use a sine flow where the phase is randomized at each period, leading to a flow 
that has a positive Lyapunov exponent at all stirring amplitudes and for all initial conditions. This 
avoids the coexistence of regular and chaotic regions, which would lead to inhomeogeneous mixing 
characteristics over the domain of the problem. Because of its uniform stirring, the random-phase 
sine flow is often used as a proxy for turbulent stirring at high Prandtl number, but at much lower 
computational cost. Using this flow, we study coarsening arrest, as in previous papers [TBI [2D]- 
However, we shall also use it to study mixing due to the stirring and the hyperdiffusion of the CH 
dynamics. Under this additional process the scaling hypothesis that characterizes the dynamics of 
bubble formation, namely that there is one dominant lengthscale [H [20] , breaks down. 

Having found a mechanism that counteracts domain formation and promotes mixing, we introduce 
a quantitative measure of homogenization, studying the p^^ power-mean fluctuation of the concen- 
tration about its average value. We specialize to the symmetric mixture in which equal amounts 
of both binary fluid components are present. By fluctuations about the average value, we mean 
spatial fluctuations around the (constant) average spatial concentration, which we then average 
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over space and time. In effect, we study the time-averaged norm of the concentration. If this 
quantity is small, the average deviation of the system about the homogeneous state is small, and 
we therefore use this quantity as a proxy for the mixedness of the fluid. An approach similar to 
this has already been taken in the theory of miscible fluids [281 ISH [30l [31]. There, the equation of 
interest is the advection-diffusion equation, and fluctuations about the mean are measured by the 
variance or centred second power-mean of the fluid concentration [321 [33] . The variance is reduced 
by certain stirring mechanisms. By specifying a source term, it is possible to state the maximum 
amount by which a given flow can reduce the variance, and hence mix the fluid. By quantifying the 
variance reduction, one can classify flows according to how effective they are at mixing. Just as the 
linearity of the advection-diffusion equation suggested the variance as a natural way of measuring 
fluctuations in the concentration, the nonlinearity of the CH equation and its concomitant free 
energy (cubic and quartic in the concentration, respectively) will fix our attention on the fourth 
power-mean of the concentration fluctuations. Owing to Holder's inequality, a binary liquid that 
is well mixed in this sense will also be well mixed in the variance sense. The advantage we gain in 
considering the fourth power-mean is the derivation of explicit bounds on our chosen measure of 
mixedness, which have manifest flow- and source-dependence. 

Due to the relevance of phase-separating thin films in industrial applications [3l [3l] , many exper- 
iments and numerical simulations focus on understanding how phase separation is altered if the 
binary fluid forms a thin layer on a substrate. One potential application of phase separation in thin 
films is in self-assembly [351 [Ml [37] . Here molecules (usually residing in a thin layer) respond to an 
energy-minimization requirement by spontaneously forming large-scale structures. Equations of 
Cahn-Hilliard type have been proposed to explain the qualitative features of self-assembly [3S1 [SB] , 
and knowledge of variations in the film height could enhance these models. Therefore, we examine 
this application, and explain the main features of these studies, introducing a long-wavelength 
approximation of the coupled Navier-Stokes Cahn-Hilliard (NSCH) equations for a thin layer of 
fluid with a free surface. 

Several recent experiments have clarified the different regimes of domain growth in a binary thin 
film. Wang and Composto [39] have identified early, intermediate, and late stages of evolution. 
The early stage comprises three-dimensional domain growth, while the intermediate stage is char- 
acterized by the formation of wetting layers at the film boundaries, the thinning of the middle 
layer, and significant surface roughening. Due to the thinning of the middle layer, the sandwich- 
like structure breaks up and matter from the wetting layer flows back into the bulk. Thus, a late 
stage is reached, consisting of bubbles coated by thin wetting layers. This characterization of the 
evolution has been seen in other experiments [iQl HH [12] , although clearly a variety of behaviours 
is possible, depending on the wetting properties of the mixture. Our model captures the essential 
features of this evolution, in particular the tendency for concentration gradients to promote film 
rupture and surface roughening. 

In a series of papers. Das et al. [131 [lH [13 [IHl [IZj numerically investigate the behaviour of binary 
fluids with wetting. In [13] they study the wetting properties of a binary mixture in an ultra-thin 
film. The behaviour is different from that in bulk mixtures. In bulk mixtures where one component 
of the binary fluid is preferentially attracted to the boundary, a layer rich in this component is 
established at the boundary, followed by a depletion layer. This layered structure propagates into 
the bulk and is called a spinodal wave [HI [IHl [IS] • In ultra-thin films, where the film thickness 
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Figure 1.1: Spinodal decomposition with wetting. Subfigure (a) provides a schematic representa- 
tion of the spinodal wave at the lower horizontal wall. This structure forms when one of the binary 
fluid components wets the wall. Layers of alternating demixed regions propagate into the bulk 
until the wall effects are negligible. This effect is suppressed in thin films. Subfigure (b) shows the 
completely wet thin film, in which a single binary fluid component wets the horizontal bounding 
walls. The film is too thin to support a spinodal wave. Subfigure (c) shows the partially wet thin 
film, in which both binary fluid components wet the horizontal bounding walls, and there is no 
vertical domain structure. In this report we shall focus on the partially wet case. 



is less than a spinodal wavelength, this layered structure is suppressed. Two types of behaviour 
are then possible. In the partially wet case, both fluid components come into contact with the 
film boundaries. The ultimate state of the system is a domain-like structure extending in the 
lateral directions. The domains grow in time as t^^^, indicating Lifshitz-Slyozov diffusion [18] . 
In the other case, complete wetting, only one of the fluid components is in contact with the film 



boundaries. These different wetting regimes are portrayed schematically in Fig. LI Our focus in 
this report is on the partially wet case. 

The papers of Das et al. elucidate the role of wetting and film thickness on phase separation, 
although they do not discuss hydrodynamics or the effect of free-surface variations on domain 
formation. We therefore focus on ultra-thin films with a variable free surface, and for simplicity 
we restrict our attention to the case where both fluids experience the same interaction with the 
substrate and free surface. The model we introduce is based on the coupled NSCH equations. 
With an applied external forcing, the model highlights the effect of the dynamical backreaction of 
concentration gradients on the flow, a useful feature in applications where control of phase sepa- 
ration is required [37|. When the binary fluid forms a thin film on a substrate, a long- wavelength 
approximation simplifies the NSCH equations, which reduce to a pair of coupled evolution equa- 
tions for the free surface and concentration. If h {x, t) is the scaled free-surface height, and c {x, t) 
is the binary fluid concentration, then the dimensionless equations take the form 

dh d 

Vx-J = 0, -7r{ch) + V±- (Jc) = V±- {hV±fi), (1.1a) 



dt ^ ' dt 



where 



fi = c^ -c - -^V±-{hV±c). (1.1c) 
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The symbol V_l = {dx,dy) denotes the gradient operator in the lateral directions. The positive 
constants C, r, and measure surface tension, coupling strength between concentration and 
free-surface variations, and the scaled interfacial thickness, respectively. Additionally, F is the 
dimensionless, spatially-varying surface tension, and is the body-force potential acting on the 
film. In this report we take = —Ah~^, A > 0, the repulsive Van der Waals potential [IH]. This 
choice stabilizes the film and prevents rupture. Although rupture is in itself an important feature 
in thin-film equations [50l I5T1 [52] , in this report we are interested in late-time phase separation 
and it is therefore undesirable. 

The analysis of thin-film equations was given great impetus by Bernis and Friedman in |53] . They 
focus on the basic thin-film equation 

' ^h-m . (1.2) 



dt dx \ dx^ ^ 

with no-fiux boundary conditions on a line segment, and smooth nonnegative initial conditions. 
For n = 1 this equation describes a thin bridge between two masses of fluid in a Hele-Shaw cell, for 
n < 3 it is used in slip models as h , while for = 3 it gives the evolution of the free surface 
of a thin film experiencing capillary forces ^U\. Using a decaying free-energy functional, they prove 



the existence of nonnegative solutions to Eq. (1.2) for n > 1, while for n > 4, the solution is unique, 
strictly positive, and is almost always bounded in the H^'^ norm. This paper has inspired other 



work on the subject [511 IS21 ES], in which the effect of a Van der Waals term on Eq. (1.2) is 
investigated. These works provide results concerning regularity, long-time behaviour, and film 
rupture, in the presence of an attractive Van der Waals force. More relevant to the present work 
is the paper by Wieland and Garcke [56], in which a pair of partial differential equations describes 
the coupled evolution of free-surface variations, and surfactant concentration. The authors derive 
the relevant equations using the long-wavelength theory, obtain a decaying energy functional, and 
prove results concerning the existence and nonnegativity {h > 0) of solutions. We shall take a 
similar approach in this report. 

We have noted that equations of Cahn-Hilliard type have been used to model the self-assembly 
of molecules. In Ch. [7] we investigate the aggregation of magnetic particles using the Holm- 
Putkaradze model of aggregation. We compare and contrast the Holm-Putkaradze and the Cahn- 
Hilliard models. The formulation of these models through a gradient flow serves as a bridge 
between our study and more recent work in nanoscale modelling. 

The modelling of nanoscale physics has become important in recent years, both because of indus- 
trial applications [53, |5Hl EHl [60] , and because of the development of experiments that probe these 
small scales [HI] . In a series of papers. Holm, Putkaradze and Tronci [621 [631 EH ESI [661 EH] have 
focussed on the derivation of aggregation equations that possess emergent singular solutions. Con- 
tinuum aggregation equations have been used to model gravitational collapse and the subsequent 
emergence of stars [HZ] , the localization of biological populations [HHl ISSl [ZD] , and the self-assembly 
of nanoparticles [33]. These are complexes of atoms or molecules that form mesoscale structures 
with particle-like behaviour. The utility of the Holm-Putkaradze model lies in its emphasis on 
nonlocal physics, and the emergence of singular solutions from smooth initial data. Because of the 
singular (delta-function) behaviour of the model, it is an appropriate way to describe the universal 
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phenomena of aggregation and the subsequent formation of particle-hke structures. Indeed in this 
framework, it is possible to prescribe the dynamics of the particle-hke structures after collapse. 
Thus, the model provides a description of directed self-assembly in nanoscale physics [351 136] . 
wherein microscopic particles collapse under the potential they exert on each other, and form 
mesoscopic structures that in turn behave like particles. Our work in this area focusses on the 
vector-scalar Holm-Putkaradze model which, although derived elsewhere p2], has not yet been 
studied in any detail. This model provides an appropriate description of self-assembly of magnetic 
particles. 

Although the Cahn~Hilliard theory has a venerable history, studies of phase separation with advec- 
tion are more recent. Our result concerning homogenization by a chaotic flow is new. The applica- 
tion of (standard) functional analysis techniques to the advective Cahn-Hilliard equation, and the 
derivation of lower bounds on the amount of homogenization possible is also new. Although long- 
wavelength theory is common in applications, the thin-film Stokes Cahn-Hilliard (SCH) equations 
are entirely new, as are the numerical and analytical results concerning this model. In particular, 
we have identified a novel way of controlling phase separation by a combination of applied surface 
stresses and backreaction tuning. The vector Holm-Putkaradze model is a recent addition to the 
family of aggregation models and, while the equations were formulated elsewhere, our analysis of 
the coupled vector-scalar model is new. 

This report is organized as follows. In Ch. |2]we derive the Cahn-Hilliard and Navier-Stokes Cahn- 
Hilliard equations from a combination of intuition and variational techniques. We discuss the 
elements of functional analysis that we need in this report, and prove an existence and uniqueness 
result for the CH equation with passive advection. In Ch. [3] we outline a numerical study of the CH 
equation stirred by a passive, chaotic flow. We characterize the mixing process by the lengthscales 
involved in the problem, by the concentration probability distribution function, and by a simple 
one-dimensional model. In Ch. |4| we use the functional analysis results of Ch. [2| and obtain upper 
and lower bounds on the level of mixedness achievable by any flow. In Ch. [5] we derive the thin-film 
Stokes Cahn-Hilliard equations, highlighting the importance of the dynamical coupling between the 
concentration and fluid velocity. We obtain existence and regularity results for the model equations. 
Although such analysis gives very precise results, it is unable to tell us much about the qualitative 
features of the solutions. Therefore, in Ch. [6] we carry out numerical simulations of the SCH 
equations, and compare the results with experimental papers. The simple model introduced here 
explains the qualitative features observed in real phase-separating fluids. In Ch. [7]we focus on the 
Holm-Putkaradze model of aggregation, which serves as a bridge between the Cahn-Hilliard model 
we have studied, and more recent work in nanoscale physics. We examine the equations introduced 
by Holm, Putkaradze, and Tronci for the aggregation of oriented particles [621 |63l EH [65]. We 
investigate these equations numerically and study their evolution and aggregation properties. One 
aspect of nonlocal problems, already mentioned in [3H], is the effect of competition between the 
lengthscales of nonlocality on the system evolution. We highlight this effect with a linear stability 
analysis of the full density-magnetization equations. Finally, in Ch. |8| we present our conclusions. 
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Publications 

Many of the results presented in this thesis have aheady been pubhshed. We summarize this work 
in what follows. Unless otherwise stated, the principal author is L. O Naraigh and the co-author 
is J.-L. Thiffeault. 

• The numerical study of phase separation and passive chaotic advection is presented in [7T] . 

• The question, "given a passive flow in a binary liquid, by what amount can the fluid be 
homogenized?" is addressed in ^72|. There, the authors use functional analytic and numerical 
techniques to answer this question. 

• Long- wavelength theory is used to derive a set of thin-film Stokes Cahn-Hilliard (SCH) 
equations in f73j. The authors highlight the importance of the backreaction when the control 
of phase separation is required. 

• Using a combination of numerical analysis and simple dynamical systems arguments, the 
vector Holm-Putkaradze model is investigated in [74j. The first author is D. D. Holm and 
the co-authors are L. O Naraigh and C. Tronci. 

Statement of originality 

The work presented in this thesis is my own, unless indication is given to the contrary. 



Chapter 2 



General Formulation 



2.1 Overview 



In this chapter we outhne the thermodynamical formahsm that gives rise to the Cahn-Hilhard 
equation and demonstrate its couphng to a flow field. We introduce the Navier-Stokes Cahn- 
Hilliard equations and obtain evolution equations for the physical quantities associated with the 
dynamical variables. We also discuss the mathematical formalism necessary to this report. 



2.2 The Cahn— Hilliard equation 

In this report, the Cahn-Hilliard model is viewed as an order-parameter equation, derived from 
thermodynamical considerations. An order-parameter description is independent of the detailed 
nature of the system being considered [7S], and thus the equation is completely general, describing 
any two-component system where the mixed state is energetically unfavorable, and where the total 
amount of matter is conserved. Let c{x,t) be the order parameter, or simply put, the concentration 
of the binary fiuid. We write down a free energy for the system in n dimensions, as a functional 
of the concentration c and temperature T, 



F[c,T] =60 dJ'x 



2 



1 - 



T 



+ h|Vc| 



(2.1) 



where Tc is a critical temperature, and where the gradient energy density |Vc| penalizes the 
formation of sharp interfaces. Indeed, we shall see that the lengthscale ^7 represents the typical 
thickness of a smoothly-varying transition zone between unmixed regions, and that the finite 
thickness of these zones prevents the formation of infinite gradients in the problem. The constant 
Sq appears in order to give Eq. (2.1 ) the correct units: Eq = kBTc/\fl\ has units of energy density. We 



shall, however, omit this constant if possible. The eager reader will note that Eq. (2.1) is nothing 
other than the Ginzburg-Landau free energy [75]. The system admits a critical temperature Tc; 
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this is a threshold temperature across which the system exhibits quahtatively different behaviour. 
For T > Tc, the energy density 



/o(c,T) = |c^-i(l--)c^ 



T 

c 

has a single minimum at c = 0, and thus the mixed state is energetically favourable. On the 
other hand, for T < Tc, the energy density has a double-well structure, with two stable minima 



at c = ±a/1 — (T/Tc). These minima describe the pure phases of the mixture. The mixed state 
c = becomes a local maximum and is unstable (as the linear stability analysis in this section 
will demonstrate). We shall focus on the subcritical case, where spontaneous demixing is possible. 
Indeed, in the following derivation, we shall set T = 0, eliminating all thermal noise from the 
system. This simplification will be justified at the end of the derivation. 

To minimize the free energy F [c] = F [c, 0] , the concentration evolves according to an equation of 
the type 

dc _ ^6F 

where O is an operator to be determined. Moreover, we are interested in conserved dynamics; 
that is, the total concentration d'^xc [x, t) must be constant. If, as in many applications |2j, the 
system is invariant under changes of parity, the simplest form for the operator O is, 

6(-) = w-[0(c)v(-)], 

where D is a diffusion coefficient and (c) is a nonnegative function called the mobility. The choice 
of mobility function is determined by physical considerations; for many of our applications, we 
shall simply set = 1, while an in-depth discussion of the choice of functional form of is given 



in Sec. 3.7 Under these assumptions, the Cahn-Hilliard (CH) equation is 

'g=DA{c'-c- 7Ac) . (2.2) 

The order-parameter nature of this equation has resulted in its use in situations far beyond its 
original purpose as a model of the phase separation of a binary alloy [2J: it has found applications 
in polymer physics [23], in interfacial flows [19], and in mathematical biology [76]. By identifying 
the chemical potential [77] fi = — c — 7AC and the flux 



'CH 



-DVfi, (2.3) 



the equation (2.2) can be recast in flux-conservative form: 

dc „ ^ 

3j + V.Jo„ = 0. 

SO that, upon applying Gauss's theorem, we have, 

^ f d^xc (x, t) + f dAnfi ■ JcH = 0, 
dt Jq J an 
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where is the problem domain in n dimensions, is its boundary. The quantity dAn is the 
differential element of area and n is the outward-pointing unit normal. By specifying the boundary 
conditions (BCs) 

n • Vc = n ■ JcH = 0, on dn, (2.4) 

the total concentration is conserved, and when combined with initial data, this gives sufficient 
information to solve Eq. (2.2). 



We verify that the concentration evolution Eq. (2.2) causes the decay of free energy, as desired. 



dF 
~dt 



-D / rf"x|V/i|' + 



dAr 



an 



Oc 

Dun ■ Vu + —h ■ Vc 

at 



Using the BCs (2.4), this is simply 



dF 
~dt 



-D / d"a;|V/i|^ 



and the free energy is a decaying function of time. 

Let us investigate the linear stability of a constant state c (a;, t) 
analysis of Eq. (2.2) about this state, we obtain an equation for the perturbation (5c, 

d 



(2.5) 



Cq. Performing a linear stability 



dt 



5c = {3cl - 1) DA5c - -fDA'^Sc, 



and by making the ansatz 6c{x,t) oc e^'^'^e""*, we obtain an expression for the growth rate of 
perturbations. 



a 



-De[{3ci-i)+^e]. 

For —l/^/S < Co < perturbations on scales below the critical lengthscale 



A 



crit 



1 

2^ 



7 



\3cl 
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decay, while perturbations on scales above Acrit grow. If the initial data contain a spectrum of 
lengthscales, this stability analysis demonstrates that small scales decay in amplitude, while large 
scales grow, suggesting the formation of large-scale structures. In later chapters we shall see that 
this is indeed the case. This linear stability analysis yields another piece of information: if the box 
size L of the problem is less than the critical scale Acrit, then no large-scale structures can form, 
and the perturbations always decay to zero. This is the statement that the problem dimensions 
are so small as to favour the hyperdiffusion of Eq. (2.2) over phase separation, a useful result to 
keep in mind when considering phase separating liquids in ultrathin films. 



The existence of equilibrium solutions to the CH equation (2.2) is well known [7S], and we provide 
the details here because they will be useful in the following chapters. It is known that the late-time 
solution is one of bubbles, for which c ~ ±1, with transition regions of thickness y/j separating 
these bubbles. It is possible to derive an equilibrium solution that describes these transition zones. 
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We write down the equilibrium version of the Cahn-Hilhard equation, 

- c - 7AC = 0. 

Let us set up a coordinate system on the bubble boundary with coordinates both perpendicular 
and parallel to the boundary. If the bubble's radius of curvature is large, as indeed it will be at 
late times, then the equilibrium CH equation reduces to 



c — c — 7 



0, 



(2.6) 



where g is the coordinate perpendicular to the bubble interface (we discuss the effects of a finite 



radius of curvature in Sec. 3.6). The boundary conditions place a single isolated bubble (c = +1) 
in a sea of negative concentration, c = — 1. Thus, the boundary conditions aX g = ±00 are c = ±1, 



c' = 0. Multiplying Eq. (2.6) by the integrating factor dc/dg and using the boundary conditions, 
we obtain the equation 

7 / dcV 
2 \dgj 

with solution 

(2.8) 



/o(c) 



(2.7) 



c (g) = tanh (^g/ a/27 j 



consistent with the given boundary conditions. Note that for 7-^0, Eq. (2.8) becomes a step 



function, with delta-function first derivative. Such infinite gradients are undesirable in our work, 
and we shall therefore always take 7 > 0. 

In n dimensions, the surface tension is the free energy associated with the bubble interface, per 
unit interfacial area. That is, 



£0 



dg 



/o(c) 



7 
2 



dg) 



Upon substitution of the solution (2.8), the surface tension reduces to T 
[r] = [Energy] and F has the correct dimensions of surface tension. 



£0^8779. Thus, 



Finally, since g is the coordinate perpendicular to the interface, the gradient Vc has the form 



gdc/dg, where ^ is a unit vector perpendicular to the interface. Thus, it follows from Eq. (2.7) 
that 



/o(c) 



7 
2 



Vcl 



and the kinetic and potential energy densities are in equipartition near the transition region. 
Moreover, far from the transition region, c = ±1, and thus local equipartition also holds in the 
bulk. Integration over the domain of the problem therefore yields equipartition of kinetic and 
potential energy. 



We examine briefly the evolution equation at finite temperatures. At finite temperatures, we 
are forced to take into account thermal fluctuations in the concentration field. Equation (2.2) is 
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modified by the addition of a noise term, 

dc 

— = DA{c^ -c--fAc) +^{x,t), (2.9a) 

where C, {x, t) is a white-noise term satisfying the relation 

(x, t) ^ {x, t')) = -2DkBTA6 {x - x') 5 {t - t') . (2.9b) 



Equation (2.9) is the Cahn-Hilhard-Cooke equation [79]. Bray [80] shows that the scahng laws 
that characterize phase separation are the same for the Cahn-Hilliard and Cahn-Hilliard-Cooke 
equations, and that temperature is thus an irrelevant variable, at least when discussing the asymp- 
totic (i.e. as t — oo) features of phase separation. We are therefore justified in ignoring it in this 
report. 



2.3 The Navier— Stokes Cahn— Hilliard equations 

Using a combination of variational principles and physical intuition, we obtain a multicomponent 
flow model that couples fluid velocities to a concentration field. This approach is based loosely on 
the work of Lowengrub and Truskinowsky |19j . 

Consider a two-component binary liquid, with components labelled by the index i = 1,2. In a 
small volume V, the i^^ component occupies a volume V^, with Vi + V2 = V. In the same way, 
if the volume V has total mass M, with the i^^ fluid component giving a contribution Mj, then 
Ml -|- M2 = M. We introduce the real mass density pi = Mi/Vi, in contrast to the apparent mass 
density pi = Mi/V. Then, pi = Mi/V = PiVi/V = piXi, where Xi is the i^^ volume fraction. 
That is, given the characteristic function Xi {x), with Xi {x) = 1 if a; is in the i^^ phase and zero 
otherwise, we have 

X^iV) = J^ [ X,{x)rx, 2 = 1,2, 

V 



\V\ 

where X2 = 1 — Xi- With this knowledge, we introduce the binary liquid density, 

M Ml M2 
P= y = ^ + ^ = piXi + P2X2- 

To obtain a more convenient expression for p, we make use of the following expression for the 
volume fraction Xi, 

_Yi - YiEiP. - MiR. - R. 

V ~ Vp p,~ M p, 
where we have defined the mass fraction q, 

Mip 

Xipi = -jf = Cip, 1 = 1,2. 
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Thus, we obtain two equations for xi, 

piXi = cip, P = Xi (Pi - P2) + P2- 

Eliminating xi between the equations, we obtain an expression for the density p which depends 
on the mass fraction ci = c, 

i = ^ + ^. (2.10) 

P Pl P2 

This is the so-called simple-mixture equation. Note that the density so defined is equivalent to the 
identity p = Pi + P2- 

Let Vi be the velocity at which the i^^ component travels in the mixture. We require that the mass 
of each binary fluid component be conserved, 

^ + V • (vrp^) = 0. 
By identifying the mass-weighted velocity field 

^=Pl]!l±P^_ (2.11) 

Pl + P2 

we obtain a conservation law for the mixture variables {v,p = pi + P2), 

f^+V-{vp)=0. (2.12) 

In this report we focus on density-matched fluids, for which pi = P2 = constant. Then the mixture 
density is simply 

1 _ 1 _ 1 
P Pl ' Po' 



and thus the mass- conservation law (2.12) reduces to the incompressibility condition 

V-i; = 0. (2.13) 
Moreover, we postulate an evolution equation for the concentration (mass fraction) c(a;,t), inde- 



pendently of the analysis in Sec. 2.2 



p{^^+v-V^ =DApD (2.14) 

where pn is a chemical potential to be determined, and D is a constant with the dimensions 
[Length]^"" [Mass] [T]~\ 

To find the dynamics of the concentration and velocity fields, we momentarily set D to zero, and 
assume that the concentration c{x,t) is in fact conserved; that is, it satisfies 
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In the Lagrangian picture of fluid mechanics, this expression has the form 

/ a 

where the partial derivative {d/dr)^ indicates time differentiation with respect to a flxed La- 
grangian label a, and where the fluid parcel with label a has the velocity 



I =V[X,T). 
a 



We write down a Lagrangian for the fluid velocity and phase-separating concentration, 



L= \ I (Tai - 



^^""^ [ [/o (c) + i7|Vc|2] + [ rafl ^ ^ ^ 

Jn Jn \ 



KdrJ Po Jn ^ Jn \ PoJ Po 

where d"'a = dm = pd'^x, and a is a Lagrangian label. The density p is then given as the Jacobian 

(9(ai, ...,a„) 



P 



d{xi, ...,Xn)' 



and the term J^d'^a[l — (p/po)] (p/Po) is the Lagrangian multiplier that enforces the constant- 
density condition. The Lagrangian action is then 



rt2 

S= drL. 
Jti 

Following standard practice [HI], stationarity of this action with respect to variations in the La- 
grangian trajectories x (a, r) yields the equations of motion. The variations along trajectories take 
place at flxed Lagrangian label, and the variation in the density has the form 

6p = — pV • 6x, 

which follows from the constancy of pd'^x along particle trajectories, [dr {pd"'x)]^ = 0. Moreover, 
the related variations in concentration have the form, 

6c = Vc ■ Sx, diSc = (didjc) Sxj + djC (diSxj) , SdiC = {didjc) Sxj, 

and hence 

S (dic) = di {5c) — djC (diSxj) . 

Throughout the report, unless indication is given to the contrary, we use the summation convention 
and sum over repeated indices. 

Using these relations, we vary the action S over particle trajectories x (a, r), at flxed particle label. 
We choose variations Sx that vanish at the end times ti and t2, and at the container boundary 
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an. Thus, 



ss = ~ 



t2 



tl 



^d^x 1 



dt I (Ta 
'n 

£o 

Po Ju 



V 



PP 
Pi 



2\ H 



• 5x + 



[\t [d-J-lfl-i-] 



dt I d^a—di {pdicdjc) 6xj 

n P Po Jti 



t2 



dt / d^a 



f, (c) - ■ (pVc) 



5c. 



Using 6S/6p = 0, we obtain the constraint p = pQ. Then using the result 6S/6x = 0, we obtain 
the equation of motion 

= Vp - —V ■ VcVc 

Po Po 

Finally, the stationarity result 6S/6c = gives the constraint 

/ac)-7Ac = 0. 



We rewrite the equations in Eulerian form, 



^ + vWv = -— Vp- — V- (VcVc), 
dt po po 

= /i = /ac)-7Ac. 
By making the identification = P and restoring D 7^ 0, we obtain the equations 

^ + . = Vp - —V ■ (VcVc) 

dt Po po 



(2.15a) 



^ + ^.Vc=-A(/ac)-7Ac) 

Ot Po 



(2.15b) 

V-i; = 0. (2.15c) 

We observe that the Lagrange multiplier that enforces the constant density plays the role of the 
pressure in these equations. The Lagrange multiplier is not arbitrary: it is fixed by the equation 



p = -poA-' 



di {vjdjVi) + —didj {dicdjc) 
Po 



Finally, we model the effect of viscosity by adding a diffusion term to Eq. (|2.15a|), which gives rise 

(2.16a) 



to the following constant-density Navier-Stokes Cahn-Hilliard (NSCH) equations 

dv 



dt 



V ■ Vv 



1 £nT 

— Vp - —V ■ (VcVc) + uAv, 
Po Po 



^ + ^.Vc=-A(/ac)-7Ac) 
ot Po 

V ■v = 0, 



(2.16b) 
(2.16c) 
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where u is the kinematic viscosity of the mixture. By specifying the viscosity term as we have 
done, we confine ourselves to Newtonian fluids. These equations possess an energy E (t), 

E = ^ [ rxv^ + 60 [ dTx [/o (c) + ^7 |Vc|'] , 

which is dissipated according to the relation 

^ = -u [ rx \Vv\' - ^ / |V/i|' < 0. 
dt Jn po Jn 



We comment on the symmetries of Eq. (2.16). By identifying a tensor 



T. - -i^A 

1 ij 0, 



Po 



f dvi 



\ dxj dxi ^ 



£07 dc dc 
po dxi dxj ' 



Eqs. (2.16) can be rewritten as 



dv 



dc 
dt 



+ i; . VV = V ■ T, 



D 



^.Vc=-A(/ac)-7Ac) 



Po 



(2.17) 

(2.18a) 

(2.18b) 
(2.18c) 



c^ — c as in Sec. 



2.2 



V -i; = 0. 

this set of equations is manifestly invariant under the change 



Choosing f'^ (c) 

symmetry present in the basic CH equation. The form of the tensor in Eq. (|2.17l) is 



essential for preserving this symmetry. One way to picture this symmetry is to imagine that 
the concentration c(a;,t) describes a binary mixture of a red and a green fluid, with c = +1(— 1) 



indicating a region of pure red (green). The invariance of the equations (2.18) is the statement that 
a person with dichromatism can be trusted to perform experiments on the system. Furthermore, 



owing to the form of the tensor in Eq. (2.17), the equations (2.18) are rotationally invariant. 



Indeed, Eqs. (2.18) provide a simple model in which these symmetries hold, and this reasoning 



provides another route to obtaining model equations for advection and phase separation. 



We note that the passive version of Eq. (2.18) is obtained by setting the backreaction term in the 



tensor Tij to zero. In this limit, it is possible to imagine the incompressible velocity field v {x,t) 
as being prescribed, and we consider the much simpler equation 



dc 
dt 



+ vVc = DA (/o (c) - 7Ac) , V ■ = 0, 



where D = D/po is a diffusion coefficient. We shall be concerned with this passive advection 
equation in Chapters |3] and |4| 
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2.4 Inequalities: Relations between function spaces 

In this section we outline some of the basic tools necessary for the mathematical analysis of partial 
differential equations (PDEs). We shall study a periodic domain fl in M", although these results 
require only that Q be compact with Lipschitz boundary dQ. We consider classes of functions that 
form a vector space with respect to a given norm. By deriving inequalities for the norms, we shall 
obtain relations between the vector spaces of functions. 

We study functions of type 



The norm of the function / is the number 



i/p 

d-x\f\P] , p>l, (2.19) 



provided the integral exists. If the integral does exist, we write / G L^(fi). The set (Q) is 
closed under addition and scalar multiplication of functions and therefore forms a vector space. 
This closure result will be useful in Ch. [s] We also have the Sobolev {q,p) norm of the function /, 

i/p 

d-x T iD^f 1 , g > 0, p > 1 (2.20) 



n 



n+---+o!n<q / 



provided the integrals exist. Here D" refers to the derivative 

^oi+...+a„ 



dx?\..dx: 



n 



If the integrals in Eq. (2.20) exist we write / G if '^(f2). Again, it can be shown that if^'^ (fi) 



forms a vector space, called the Sobolev (g,p) space. Finally, for functions of the type 

/ : X [0, T] ^ M, 

we introduce the norm 

provided the integral exists. If the integral exists we write / G U' {0,T; H'^'^ (Q)), and identify 
another function space, L'" (0, T; iJ'^'^ (^))- 

To make progress in deriving relations between these norms, we introduce the following inequalities 
for real numbers. 

• Triangle inequality. If a and b are real numbers, then 

\a + b\ < \a\ + \b\. 
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• Young's first inequality. If a and b are real numbers, and if k, is any positive number, 
then 

2 

\ab\ < Ka + —. 
Proof: Since the square of any real number is nonnegative, 

6^ 



o<(v^H-^|6|)' 



Ko^ + \ab\, 



and the result follows. 



Young's second inequality. If a and b are positve real numbers, and if p, g > 1 satisfy 
1/p + 1/q — 1, then 

, aP h" 

ab< \ . 

P Q 

Proof: Using the string of identities 

flog(aP) log (6*1 

ab — exp (log {ab)) — exp (log (a) + log (6)) = exp — - H — - 



P 

and the convexity of the exponential function, e*^+(^~*)2/ < _|_ _ ^y^ ^qj. ^ ^ 
obtain, by identifying t — 1/p, 1/q — 1 — t, 

, aP bi 

ab< h — , 

P Q 

as required. 

Using this knowledge, we prove some results for function spaces. 

• Holder's inequality. Let / G Lp (il) and g G L'' {^), with p,q > 1 and 1/p + l/q = 1. 
Then 

mil < wfWpMU- 

Proof: Without loss of generahty, we shall prove the inequality for the case when = 
llg'llg = 1. Then it suffices to show that < 1. By Young's second inequality, 

1/91 < ^ + ^. 

p q 

Integration over the domain fl gives the required result. 

• Cauchy-SchwEirz inequality. When p = g = 2 in Holder's inequality, we have the result 

WfaWi < WfhWgh. 
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Holder's inequality for p = oo, q = 1. If the function / is in tlie class L°° (Q), that is, if 

ll/lloo = ess-sup^.gf^|/ (x) I < oo, 
then Holder's inequality holds with p = oo and q = 1, 

\\f9\\i<\\f\\oo\\9\\i. 

Monotonicity of norms. By setting g (x) = 1 in Holder's inequality, and / = l^l*", where 
r > 1, we obtain the result 

< ll^ll ^<^_ 



A function that is bounded in the [Q) sense is therefore bounded in the L'^ {fl) sense. Such 
a result is called a continuous embedding |H2]- We write, 

L' {Q) C {Q) r<s. 

Minkowski's inequality. Let f,g ^ L^(f^). Then, 

\\f + 9\\p<\\f\\p+\\9\\p. 
Proof: Using the triangle inequality, we have. 



rx\f + gf< / rx\f + gry\+ / rx\f + gr'\g\ 
n Jn Jn 

Applying Holder's inequality to this relation, we have, 
rx\f + g\P 



a, 



< / d^x\f + g\' 



n 



d-x\f\^ ] + / d-x\g 



n J \JQ 



IP 



and the result follows. 

Poincare's inequality. Let the function / be a smooth, square-integrable function on the 
periodic domain Q = [0, L]". Then we have the relation 



ll/-/ll.<(^j liv/ll^. 

where / = |^^|^^ d"'xf (x) is the mean value of the function /. 

Proof: Let m be a variable running over all values in Ng = (N U {0})". By Parseval's 
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identity, 



27r 



2 



/ 



Gagliardo— Nirenberg— Sobolev inequality. Let the function / be a mean-zero, smooth, 
square-integrable function on the periodic domain VL = [0, L]^. Then for 1 < q,r < oo and j 
and m integers satisfying < j < m, there is a constant < ci < oo such that 

\\D^f\\p<ci\\n^fUf\\'r, 

where 

1 j ( \ m\ 1 — a 

- = — ha 



p n \r n J g ' 



and for a in the interval 



< a < 1, 
m 



with the restriction that if m — j — n/r is a nonnegative integer, then a must equal j /m. 
Proof: The proof is given in |83j . 

Gagliardo— Nirenberg—Sobolev inequality (celebrated special case). Let the function 
/ be a smooth, square-integrable function on the periodic domain f2 = [0, L]". For j = 0, 

m = 1, a = ^, r < n, and 

p = q = = r > r, 

n — r 

the Gagliardo-Nirenberg-Sobolev inequality reduces to 

||/-7||r* <Ci||V/||., 

where r* is called the Sobolev conjugate of r. Rewriting this 

\\f\\Lr'in)<C2{\\Vfhr^n) + \f\), 

we have the embedding 

H^''' (fi) C U'* (n) . 
Here Ci and C2 are constants that are independent of the function /. 



Using the Minkowski inequality, and the result 

= iff / = almost everywhere. 
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it follows that || ■ ||lp(q) and || ■ || (n) are norms on the vector spaces (Q) and H'^'^ (Q) respec- 
tively. These normed vector spaces are, moreover, complete, and therefore form Banach spaces [8^ . 

As a conclusion to this theory section, we provide a list nonstandard inequalities that will be vital 
in Ch. [5] in analyzing the thin-film Stokes Cahn-Hilliard equations in one spatial dimension. 



Let Q he & periodic domain in M and let f -.n-^R belong to the class H'^'^ (n). Then the 
following inequality holds, 

sup I/I < l||/||i + ||/,||i. (2.21) 
n 1^ 

Proof: Using the Fundamental Theorem of Calculus, we have 

Let 1/ (x) I be the maximum of the function |/| over the set Vt. That is, the maximum value 
of I /I is attained at x. Then 



sup I/I < |/(a)| + / ds 



df 



ds 



< |/(a)|+ / ds 



df 



ds 



Since this is true for any a G f2, by integrating over a, we obtain the inequality 

1 



sup I/I < 

< 



1 



2 + VLll/, 



lx\\2, 



as required. 

Let / : ^ M belong to the class H'^'^ {fl). Then the following inequality holds, 

4 



ll/.||^<^ll/..||f + Tll/llill/= 



xa; 1 • 



(2.22) 



Proof: We have the identity J^f^dx = —J^ffxT, true for any function / with periodic 

boundary conditions (in general, this identity is true when ffxlan = 0). With Holder's 
inequality, this yields 

||/.||^< II/II00II/..II1. 



Using the relation (2.21), this becomes 

ml < 



< 



i + ii/xiii 

L\\f.\\2 



1 1 fxx 111; 

1 1 fxx 111) 



2.5. Galerkin approximations 



37 



which is a quadratic inequahty in ||/x||2, with solution 



ll/.||2<| 



By sacrificing the sharpness of the bound, we obtain a simpler one, 



||/.||^<^||/..||? + tII/||i||/.x.||i, 



as required. 



We shall use some of this formalism in the next section in proving the existence and uniqueness of 
solutions for a given PDE. 



2.5 Galerkin approximations 



In this section we study the typical fourth-order partial differential equation 



du 
di 



A\ = f{x,t) 



(2.23) 



on the periodic domain [0, L]" in M", although the discussion is equally valid for any smooth 
bounded domain together with appropriate boundary conditions. The function / (x, t) belongs to 
the class (0,T; {^)), for any T > 0. Although the equation that forms the subject of this 



report is nonlinear, the intuition we gain in analyzing the linear equation (2.23) will enable us to 



study the advective Cahn-Hilliard equation (3.1), and to prove the uniqueness result of Ch. |5} 
We specify the initial condition 



u{x,0) =Uo (x) e H'^^ {VL) 



and introduce the pairing 



IM, fj= / (r'Xu{x) V {x) . 



(2.24) 



The strategy we take in analyzing Eq. (2.23) is to define weak solutions, construct a sequence of 



approximate weak solutions, show that this sequence converges to a limit that we identify as the 
solution, and obtain uniqueness results for this construction. 



A function u{x,t) satisfies Eq. (2.23) weakly on f2 x (0,T] if the following integral identity 
holds, 

j^{u,v) + {u,Q*v) = {f,v), (2.25) 
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V (x) is a smooth periodic test function, and Q* = Q = A^. 



We construct approximate solutions of Eq. (2.23) as Galerkin sums, 



N 



(2.26) 



i=0 



where the 0i's are eigenf unctions of the Laplace operator on Q, with eigenvalues — A^, and 
where 0o is the constant eigenf unction. The set {0o, 0Af, •••} is a complete orthonormal 
basis with respect to the pairing (2.24). If the equation (2.23) has initial data 



u {x, 0) = Mo (x) = ^ 0i (x) 7]] 



then the approximate solution m^t {x,t) is seeded with the initial data 



N 



UN (x, 0) = u% (x) = ^ 0i (x) 



1=0 



and u% (x) converges strongly to Uq (x) in the norm. 



To obtain the ?7j's, we demand that the Galerkin approximation (2.26) satisfy the following 
weak form of the PDE, 

d 



dt 



(mat, 0i) + (mat, Q*(t)i) = (/, 0i 



(2.27) 



for (pi G {(pQ, (piy}. This is simply 

drji 
dt 



N 
j=0 



31 Yi) 1 



where the /j's are given by the relation / = Yl'iLo fi (^) 0* {^)- matrix form, this is 

drj 



dt 



f-Ari, 



(2.28) 



where Aij = {Q(f)i, In view of the boundedness of the operator Q, the right-hand side 
of this expression is globally Lipschitz in the vector 77, and thus local existence theory [85] 
guarantees that Eq. ( |2.28 ) has a unique solution that is Lipschitz continuous in time, in a 
small interval [0,0"). Thus, we have prescribed un {x,t) in a small time-interval [0,cr). 

To continue this solution to another time interval [a, cxi), and thus to the whole interval 
[0,T), we must find a bound on ||MAr||2 that depends only on the initial data. To do this, let 
us first of all replace (pi with un in Eq. (2.27), and obtain 

id 



'dt 



\\un\\1 + \\Aun\\1 < 



1 
4k 



2 , .11 ||2 
2 + II II 2) 



(2.29) 
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where k is an arbitrary positive constant that comes from Young's first inequahty. Hence, 



[1—2k,T) i sup ||MAr||2(t)) < sup ||MAr||2 (0) + 7^ / 

Yte[o,T] J Are[o,oo) 2k Jo 



2 (s) ds. 



We choose k such that 2kT < 1. Then, 

sup \\uN\\lit) < ^— — 

te[o,T] (1 - 2kT) 



1 

sup \\un\\1 (0) + — / 

^GfO,oo) Jo 



Ne[0,oo) 

<t<a <T, 



2 (s) ds 



< ki < oo, 

where the number ki depends only on the initial conditions, and on the time T. 



(2.30) 



Additional bounds are obtained by replacing 0j with A^w^r in Eq. (2.27). Then we obtain 
the inequality 



^ d 



\\Aun\\1 + \\A\n\\1 < 



1 



+ /tllA^-Uivl 



2? 



(2.31) 



where k is an arbitrary positive constant that comes from Young's first inequality. Choosing 
K, = 1 gives 



sup \\AtiN\\lit) < sup \\Aun\\1{0) + 1 



te[o,o-] 



< k2 < oo, 



2 (s) ds, 



< a < T, 



(2.32) 



where the number /c2 depends only on the initial conditions, and on the time T. Choosing 
K = I gives 



[ \\A\N\\l{s)ds< sup ||Aujv||2(0) + / 

Jo Afe[0,oo) Jo 



2 (s) ds, 



0<t<a <T. 



(2.33) 



The uniform bounds (2.30), (2.32), and (2.33) enable us to continue the approximate solution 



UN{x,t) to the whole time interval (0,T]. Moreover, a uniformly bounded (in (f2)) se- 
quence of functions {un}n=o contains a subsequence (still denoted by {un}) which converges 
weakly to a limit u, in the following sens^ 



as ^ oo, V E L'^ (f2) . 



Since the bound (2.30) is independent of N and a, it holds in the limit mat 



\\u\\Ut) < k^ < oo. 



<t<T. 



*A bounded sequence on a Hilbert space always has a weakly convergent subsequence. Moreover, for weak 
convergence on a Hilbert space, Xn ^ x, we have ||a;|| < liniinf„_^oo ll^Jnll- 
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Thus, by Poincare's inequality, we have a hierarchy of weak convergence results. 



veL^ in) 



dt {\/Au]\i,ip) dt (Vm, v^) 



[ dt (A\Ar,v?) ^ f dt (Am,v?) 
Jo Jo 



as 



A^->oo, ^ e L^O,T;L'^ (Q)) . 



Since the bounds we have obtained are independent of A^, we have the conditions V Am, A^u G 
(0, T; (fi)), in addition to the fact that u, Vu, Au e (0, T; (fi)). 



By multiplying the Galerkin weak solution 

d 



dt 



UN, + (un, Q*(j)i) = (/, (f)i) 



(2.34) 



by a function r (t) that vanishes at t = 0, T and integrating the result over time, using the 
limiting relation un ^ u, we obtain the relation 



dt (m, <Pi)-77+ / dt (m, Q*<Pi) T{t)= / dt if, 0i) r (t) , 
^'^ Jo Jo 



(2.35) 



By superposition, we can replace (pi by a smooth function v (x). Then, taking r (t) to be a 
distribution on (0,T), we obtain the result 



d 
dt 



{u,v) + {u,Q''v) = {f,v) 



(2.36) 



which is Eq. (2.25) 



Inspection of Eq. (2.25) shows that {u,v) is in the class H^'^ ([0,T]), for any smooth v (x) 
and it is therefore continuous. It follows that u {x, 0) = Uq (x), as required. 



Extra regularity is obtained by replacing the function v in Eq. (2.27) with du^/dt. Then 



T 



dt 



dt 



< sup \\Auj,\\l{0)+ / dtWfWl 

2 Ne[0,oo) Jo 



To characterize the solution further, we make use of the theory of semigroups [86]. Formally, a 
solution to the equation 



du 
di 



+ A^u = f {x,t) , u (x, 0) = uq (x) 
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can be written as 



The object Pt := e 
semigroup property 



e uo 



Jo 



(2.37) 



tA-" 



t > is a map from the Banach space {Q) to itself, and possesses the 



PtPs = Pt+s, t,S>0. 



If Mo G H'^'^ (Q) and if / (■, ■) is C°° in its argument^ then it can be shown [HSl EZ] that the 
solution (2.37) is a unique, strong solution of the equation (2.23), and coincides with the solution 
constructed by the Galerkin approximation. Indeed, the semigroup equation (2.37) holds for 
/ = (x, t, Vu, Am, VAu), where /(■,■,■,■,■) is a function of its arguments. In this manner, 
we have obtained the unique, strong solution to the equation (2.23), which satisfies the prescribed 
initial conditions and belongs to the class 

u G (0, T; H^'^ (fi)) n (O, T; H^'^ (fi)) n H^'^ (O, T; (fi)) . 

We shall use this result, and this method of constructing solutions, throughout this report. 



2.6 A strong existence theory for the advective Cahn 
Hilliard equation 



In this section we shall prove that solutions to the equation 

dc 



dt 



+ vVc = DA (. 



c — c 



7 Ac) 



exist and are unique, once a sufficiently regular starting condition and velocity field v {x,t) are 
prescribed. We work on a periodic domain f2 = [0, L]" in M", although this restriction can be lifted. 



We shall use the machinery developed in Secs. |2.4| and |2.5[ We present this analysis for two reasons. 
Firstly, it is important to develop understanding of the basic advective Cahn-Hilliard equation 
before carrying out further numerical and analytical studies. Secondly, since the results of this 
section are a combination of old analysis of a new equation, this section serves as a bridge between 
the well-known material of this chapter, and the new results that follow. To our knowledge, no 



analysis of the advective Cahn-Hilliard equation (2.38) exists in the literature, although the proof 



presented here is similar to Elliott's and Zheng's [6J for the nonadvective case. The introduction of 



advection necessitates a new ingredient in the proof, which we provide in (2.42). We shall provide 



analysis of the long-time behaviour of Eq. (2.38) in Ch. 4 Here prove the following statement. 



weaker condition on the function / (•, •) will suffice: if / (•, •) is in space and Holder continuous in time, 
then Eq. (2.371 is still the unique, strong solution of the equation (2.23). Relaxing the assumption on / even further 



gives rise to a less smooth solution. 
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Given the advective Cahn-Hilliard equation 
dc 

— + vVc = DA{c^-c- 7Ac) , V ■ v = 0, (2.38) 

defined on the periodic domain Q = [0, L]", where D and 7 are positive constants, where 
V {x,t) is bounded in space and time in the following sense, 

and where the initial data are defined by the identity 

c (x, 0) = Co (x) G H^'^ in) , 



there exists a unique strong solution to Eq. (2.38) that lives in the class 
c G (0, T; /f2,2 p ^2 ^4,2 p ^1,2 ^2 

/or any T > 0. 



(c, m) + (c, 7L)A2n - V ■ Vn) = D (c^ - c, Au) , (2.39) 



We begin by constructing a Galerkin approximation to the weak problem 

d 
dt 

where u [x) is a smooth function on Vt. We form the Galerkin approximation 

N 
i=l 

with initial data 

00 N 

Co (x, 0) = ^ 7]^(j)i (x) , ctv (x, 0) = c^ (x) = r]°(f)i (x) , 

i=0 i=Q 

where the 0i's are the eigenfunctions of A on the domain Q, with corresponding eigenvalues — A^; 
00 is the constant eigenfunction. The Galerkin approximation cn (x,t) is then required to satisfy 
the equation 

(cat, (f)i) + {cN,'jDA'^(f)i - V ■ = D{c\- Cat, A0i) , 



dt 

or in vector form, 



= 'Y^ijVj - DX^ VjVkViBijki, (2.40) 

j=0 j,k,l=0 
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where Ay is the matrix 

Aij = {(f)i,-DA(f)j+-fDA^(f)j-v-V(f)j) , 
= {(pi, D\](l)j + '^DX^-cpj - \/^v ■ kj(t>j) , 
= Sij {D\^ + — ■ fcj) , no sum over j, 

and thus \\A\\2 exists (of course, any hnear operator on a finite-dimensional vector space has finite 
norm). On the other hand, we have the totally symmetric tensor 

Bijki = / (fxcpiCpjCpkCpi, 
Jn 

which is bounded in sense that supj^-;.; \Bijki\ < oo. It follows that the function J2jki^ijkiViVjVi is 
Lipschitz in the vector {rjo, ...,r]i\f), and local existence theory then guarantees that a solution to 



Eq. (2.40) exists in a small time-interval (0,cr) 



Next, we obtain a uniform bound on {x, t). Testing the equation (2.39) with v = cn, we obtain. 



^^\\cn\\1 + 7D\\Acn\\1 = D{c%-cn,Acn) , 

= {Vcn,Vcn)-{3cI,\Vcn\^) 

< D\\Vcn\\1 

< D\\AcNh\\cNh, 

< 7^l|Ac^ll2 + ^l|c^ll2, 
where the last line makes use of Young's first inequality with k = 'j. Thus, 



|c^ll^(T) < I sup ||c^||^(0)le^^/2^ 

Afe[0,oo) 



r dt\\Ac^\\l < ( sup llc^ll^(O) I (2.41) 
Jo yAfe[o,oo) J 

where < T < a. 

Now we consider the Galerkin approximation 

^ (cat, (f)i) + {cN^iDA^cf), - V ■ = D {c%- cn, Acpi) , 

and the chemical potential fiN = c% — cn — 'jAcn, and replace the test function 0j by fiN- Thus, 

'dcN 



-Q^^I^N ] +{v ■ Vcat, /iAf) = -D (V/iAT, VfiN) 
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or 

+ {v ■ Vcn, Hn) = -D {ViiN, V//iv) , 

where = jQd'^x ^ (c% — 1)^ + ^7 |Vc7v|^ is the approximate free energy. We examine the 

term {v ■ Vcn,IJ^n)- The parts of /in that are powers of cn drop out of this expression, owing to 
the antisymmetry of the operator v • V. We are left with 



{v ■ Vcjv, Acjv) = / (didiCN) (vjdjCN) d^x, 
Jn 

(diCN) [di (vjdjCN)] dP-x, 



= - / {diCN) (diVj) {djCN) d'^x - / (9iCjv) {v ■ V) {diCN) d'^x, 
Jo, Jn 

— — wWw'^d'^x, w — Vcat, Wij = \ {diVj + djVi) . 
Jn 



The quadratic form wWw^ satisfies < nmaxj^ iWi^l Hiflli, which gives rise to the in- 
equahty 

[ AcnV ■ VcNd^x < n (sup \Wij\] I |Vcjv|^d"x. (2.42) 

Jn J Jn 

Hence, 

dP 

^ + \\Vi^N\\l<nW^\\Vc^\\l, 
where Woo — ^^Pte[o,oo),n,i,j Integration over time yields 



/ dt\\ViiN\\l< sup FN{0)+nWoo [ dt\\VcN\\l, 

Jo Ne[0,oo) Jo 



and since cn e L'=° (0; T; (Q)) and Acjv G (0, T; (Q)), it follows that 

VAiiv,VAcAreL2 (0,T;L2 (f])) , 
independently of the order of the Galerkin approximation N . 
The same integration yields 

supFjv(t)+/ dt\\V^iN\\l< sup FNiO)+nWoo [ dt\\VcN\\l 
te[o,T] Jo Are[o,oo) Jo 

Since Jq^ <^^|| VctvH^ < OO; independently of A^, we have the bound 

I7SUP ||VcAr||2 < 00, 

[0,T] 

independently of A^. 
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Hence, from the theory developed in Sec. 2.5, it follows that 



• (cat, Vcat, Acn, V/iAf) ^ (c, Vc, Ac, /i, V/i), weakly in (0, T; (f^)); 

• ce L°°(0,r;i7i'2(fi)). 

We shall now prove that A^cat G (0, T; and hence that A (c^ — cat) = A/xo.at ^ 

(0,T;L^ (^))- We shall work in dimension n = 1, although the generalization to dimensions 



n = 2 and = 3 is straightforward. We multiply the advective Cahn-Hilliard equation (2.38) by 
A^cat and integrate over space, to obtain 

^IIAcivll^ + D^WA^cnWI 

= D {Aflo^N, A^Cat) - {v ■ VcAf, A^Cat) , /io,7V = n - Co^N 

< D (A/io,jv, A^cn) + Vc^lUI A^c^vlh, " + " = i 

p q 

< D||A/io,jv||2||A2c7v||2 + lli^llpllVcjvll'IIA^CTvlh, 

= D\\ [{3C% - l) Acat + 6cn |VcAr|^] ||2||A^CAr||2 + ||f ||p||VcAr||g||A^CAr||2. 

Let us consider the terms on the right-hand side. The first term is 

Z^ll (34-1) Ac^llsllA'cTvlb < ||34- llloollAcArllallA^c^lls, 

< I sup \\3c%-l\\l] f-^IIAcArll^ + ztillA^c^ll^') , 
Ve[o,T] y V4/ti / 

where < ki < oo is an arbitrary positive constant. The second term is 

6D||cAr IVcArHbll ||A^CAr||2 < 6L) || Cat || oo || Vcat || 4 || || A^Cat || 2, 

< 6DL=^/^||c;v||oo||Ac^||2||A2c;v||2, n = l, 

< 6DL3/4||c,v||oo (||c,v||2||A2c^||2)'^' IIA^c^lh, 



< 6DL||c,v||^lA 



Cat II 2 5 



^3 II 1 1 fi 1 1 A V n y 

< -^IWnWoo + '«2||A CAr||2, 

where we have used Young's second inequality and the Gagliardo-Nirenberg-Sobolev inequality 
in one spatial dimension. In particular, we have used the following Gagliardo-Nirenberg-Sobolev 
estimates, 

3 

IIVcatIU < -^^*||AcAr||2, n=l, 

and ^ 

||cAr||oo < ^^||cAr||2 + V^||VcAr||2 < OO, n = l] 



the procedure for obtaining these results is sketched in Sec. 2.4 The same approach, with appro- 



priate Gagliardo-Nirenberg-Sobolev estimates, gives the necessary bounds in higher dimensions. 
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Finally, the third term is 

1 



ll'i^llpl|Vciv||q||A cjvlh < — llvllpllVcjvll^ + /«4||A Cjvlla- 
Choosing ki (supjg[o,T] l|3c|f — 1||^) +fi;2+'t4 = ^jD, and integrating over [0, T] gives the inequality 

\\Acn\\1{T) + IjD [ dt\\A^CN\\l 
Jo 

< \\Acn\\1 (0) + 7^ ( sup \\3cl - [ dt\\AcN\\l 

4«1 \t(^[0,T] J Jo 

+ 5 ( sup \\cN\t] + / dt\\v\\l\\VcN\\l. 

1^2 Vtefo,Tl / 4fi;4 Jq 



A sufficient condition for the boundedness of ||Ac;v||2, and A^c^vlh obtained when p = oo 

and q — 2. We require that 



Jo 



lo<00. 



Then we have the bound 



\\Acj,\\l{T) + l^D f dt\\A^CN\\l 
Jo 

< sup ||Aca,||2(0) + -^ I sup ||34-1||^ I /" dt\\AcN\\l 

Ne[o,oo) Yte[o,T] J Jo 

+ ^ I sup llc^IlL I + ( sup llVc^lin [ dt\\v\\l,. 

1^2 \te[0,T] J "il^A \tG[0,T] J Jo 

Owing to the A^-independence of these results we obtain further information about the regularity 
of the weak solution c = limjv^oo cn, 

• AcG L°°(0,T;L2(fi)), 

• A^ce L2(0,r;L2(n)). 

It follows that A/io G L'^ (0,T;L^ (^)), which paves the way for our key result. This result also 
holds in dimensions n = 2 and n = 3. 

We obtain our final result by rewriting the advective Cahn-HiUiard equation as 

r)c 

— + -fDA^c = DA (c^ -c]-v-Vc = f. 

Ot V ^ ✓ 
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We have shown that the function / resides in (0,T;L^ (^))- Using the result of Sec. 2.5 
follows that the solution c{x,t) is in fact strong, unique, and belongs to the class 

c e (0, T; H^'^ (n)) n (O, T; H^'^ (fi)) n H^''^ (O, T; (fi)) . 



it 



Although the results in this section hold for any finite time-interval [0, T], they tell us nothing about 
the late-time behaviour of the solution to Eq. (2.38). Indeed, the exponential-in-time bounds (2.41 ) 
are distinctly unhelpful in answering this question. In Ch. |4] we shall prove the existence of long- 
time averages of the free energy F. Let us, however, take this existence on trust for now, and 
perform numerical simulations on Eq. (2.38) for the case where v {x,t) is a chaotic flow. 



Chapter 3 



Chaotic stirring of a Cahn— HiUiard fluid 



3.1 Overview 

In this chapter we study the advective Cahn-Hilhard (CH) equation 

uc 

— + vVc = DA(c^ -c--fAc) , V-v = 0, (3.1) 
ot ^ 

for the case where v {x,t) is a prescribed chaotic flow. We focus our attention on the symmetric 
mixture, in which equal amounts of both binary fluid components are present. This special case 
involves a wide range of phenomena [HI [HI 1201 EHl En] and we are therefore justified in considering 
it here. However, the results in this chapter do not depend on this assumption. By numerically 
simulating a chaotic flow by the random-phase sine map, we characterize the effects of chaotic flow 
on the late-time concentration morphology. The key tool in this study is the use of scaling laws 
derived from dimensional analysis. 



3.2 Scaling laws for the advective Cahn— Hilliard equation 

In this section we discuss the notion of dynamical equilibrium for the CH equation without flow 
in terms of the structure function. We then examine the lengthscales that arise in the presence of 
flow. 

The CH equation without flow exhibits dynamical equilibrium in the following sense. Starting 
from a concentration fleld fluctuating around the unstable equilibrium c {x, t) = 0, the late-time 
concentration fleld has properties that are time independent when lengths are measured in units of 
typical bubble size i?b- Such properties include the structure function [90] which we now introduce. 
A measure of the correlation of concentration between neighbouring points, given in Fourier space. 
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is the following: 

S(k,t) = -^ [ (Tx [ (Tx'e-'^-'^ \c {x + x') c {x) - ^'1 , (3.2) 

where (■) denotes the spatial average. We normalize this function and compute its spherical average 
to obtain the structure function 

s{k,t) = ——^— ^ (3.3) 
(27r) c2 - c 

The spherical average (k) any function (fe) is defined as 

4>{k) = — [ dcon^ik), (3.4) 

t^n J 

where dun is the element of solid angle in n dimensions and Un is its integral. Thus is the 
spherical average of the Fourier coefficient Cfc. For a symmetric binary fluid c = 0, the structure 
function is simply the spherically-averaged power spectrum: 

sik,t) = j—^^^, fore =0. (3.5) 

The dominant lengthscale is identified with the reciprocal of the most important /c-value, defined 
as the first moment of the distribution s {k,t), 

^ J^Js{k,t) dk ^ ^g^g^ 
s {k, t) dk 

Since the system we study is isotropic except on scales comparable to the size of the problem 
domain, this lengthscale is a measure of scale size in each spatial direction. That is, we have 



lost no information in performing the spherical average in Eq. (3.3). This lengthscale is in turn 
identified with the bubble size: i?b oc 1/ki. In two-dimensional simulations (n = 2) it is found 
[5l |90] that while ki and s {k, t) depend on time, k^s {k/ki,t) is a time-independent function with a 
single sharp maximum, confirming that the system is in dynamical equilibrium and that a dominant 
scale exists. Indeed, the growth law ~ t^^^ is obtained, in agreement with the evaporation- 



condensation picture of phase separation proposed by Lifshitz and Slyozov (LS exponent) 
We present some simple scaling arguments to reproduce the t^^^ scaling law and to illuminate the 



effect of flow. For v = 0, the equilibrium solution of Eq. (3.1) is c ~ ±1 in domains, with small 
transition regions of width ^ in between. Across these transition regions, it can be shown [19] 
that 

where 



r = (3.7) 
is the surface tension and k is the radius of curvature. Thus, if i?b is a typical bubble size, that is. 
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a length over which c is constant, the chemical potential associated with the bubble is 

Sfi ~ T/Rb. (3.8) 

Balancing terms in the Cahn-Hilliard equation dtc + v ■ Vc = DAfi, it follows that the time t 
required for a bubble to grow to a size _Rb is l/t ~ TD/Rl, implying the LS growth law i?b ~ 

Now the CH free energy F[c] is the system energy owing to the presence of bubbles: 

F[c]= [ rx\l{c'-lY + l^\\/cf . 

The surface tension T is the free energy per unit area. A bubble carries free energy TR^^^, where 
prefactors due to angular integration are omitted. The total free energy F [c] in the motionless 
case is then TR^~^ x A^b, where iVb is the total number of bubbles. Because the system is isotropic 
and because there is a well-defined bubble size, we estimate A^b by The free energy then 

has the scale dependence ~ ^/Rh, and using the growth law for the length i?b we obtain 

the asymptotic free energy relation 



t~^'\ (3.9) 



By introducing the tracer variance 



a2 = c2-c' (3.10) 
and restricting to a symmetric mixture c = 0, we identify cr^ ~ i7b/|^|, where fib = Ab-Rb 

n- 
b 

i?b ~ (t'^/F. (3.11) 



volume occupied by bubbles. Since F = TR? ^A^b, it follows that 



We shall verify these results in Section 3.5 Equation (3.11) will provide a useful measure of 
bubble size when a flow is imposed, since no assumption is made about the number of bubbles 
per unit volume. In a dynamical scaling regime, the lengthscales calculated from these energy 
considerations must agree with that computed from the first moment of the structure function, 
since there is only one lengthscale in such a regime 



Stirring a CH fluid introduces new lengthscales. The flow will alter the sharp power spectrum 
found above and so it may not be possible to extract definite scales from experiments or numerical 
simulations. We might expect further ambiguity of scales during a regime change (crossover), for 
example, as the bubbles are broken up and diffusion takes over. Nevertheless, for a given flow, it 
is possible to construct lengthscales from the system parameters. We shall impose a flow v {x,t) 
that is chaotic in the sense that nearby particle trajectories separate exponentially in time at a 
mean rate given by the Lyapunov exponent. 



If the advection term and the surface tension term have the same order of magnitude in a chaotic 
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flow, by balancing the terms v ■ Vc and DV^fi we obtain a scale 

R,t = {rD/\f\ 



(3.12) 



where A is the Lyapunov exponent of the flow. Because the term v ■ Vc gives an exponential ampli- 
flcation of gradients, the balance of segregation and hyperdiffusion in the term = /q (c) — 7V^c 
may be broken and hyperdiffusion may overcome segregration. This will happen on a scale -Rdiff 
given by the balance of the terms v ■ Vc and 7V^c, 



R. 



diff 



(3.13) 



In this regime, we expect mixing owing to the presence of the advection and the (hyper) diffusion. 



Using Eqs. (3.7), (3.12) and (3.13), the crossover between the bubbly and the diffusive regimes 



takes place when A ^ In the following sections we shall examine these scales and look for this 

crossover between bubbles and fllaments, the latter being characteristic of mixing by advection- 
diffusion. 



3.3 A model stirring flow 



In order to run simulations at high resolution, we make use of the sine flow model (4.25) in two 

t2 



dimensions. We are interested in the following velocity fleld, deflned on a square domain [0,L]" 
with periodic boundary conditions: 



(x, y, t) = Vq sin 
Vy (x, y, t) = Vq sin 



27ry 
L 

271X 



1 I ' 



0, 3T<t< (j + \) r, 
^0, (j + ^)r<t< (j + l)r. 



(3.14) 



where 0,- and ipj are phases that are randomized once during each flow period r and the integer j 



labels the period. We solve Eq. (3.1) with and without the flow (4.25) to investigate the scaling 



laws presented in Section 5.2 



Equation (3.1) is integrated by an operator splitting technique, in which a specialized advective 



step followed by a phase-separation step are performed at each iteration. The phase-separation 
step is the semi-implicit spectral algorithm proposed by Zhu et al. [5J for the case without flow. 
The semi-implicitness of the algorithm enables us to use a reasonably large timestep, while the 
spectral nature of the phase-separation step enables us to work at high resolution, namely 512^ 



gridpoints. By imposing the sine flow (4.25), we may use Pierrehumbert's lattice method for 



advection [271 El]- With the problem deflned on a discrete grid, this method is exact for the 
advection step. The sine flow has the added advantage that its correlation length is of order of 
the box size L, which is the largest scale in the problem. Thus, using the prediction of Lacasta et 
al. [T7] explained in Ch. [T| we expect this velocity fleld to produce coarsening arrest even at very 
small stirring amplitudes. The parameters D and 7 in Eq. (3.1) are dimensional, and we therefore 



nondimensionalize the equation before undertaking numerical simulations. Let T be a timescale 
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associated with the velocity v {x,t), let Vq be the magnitude of v {x,t), and finally, let L be the 



box size. It is then possible to write down Eq. (2.2) using a nondimensional time t' = t/T and a 
nondimensional spatial variable x' = x/L, 

^ + av-Vc = D'A' (c'-c- 7'A'c) , 
of ^ ' 

where i; is a dimensionless shape function, a = TVq/L, D' = DT/L"^, and where 7' = 7/-^^ The 
quantity D' = DT/L^ is identified with the ratio T/T^, being the ratio of the velocity timescale 
to the diffusion timescale. For ease of notation, we shall henceforth take v to mean av, and drop 
the primes from the nondimensionalized equation. Now the lattice method with its splitting of the 
advection and diffusion steps, is effective only when the spread of matter due to diffusion is much 
slower than the spread due to advection, that is, T/Td <^ 1. We therefore set D = 10~^. We shall 
also set r = 1 and L = 2tt. Following standard practice POl CS], we choose 7 ~ Ax^, the gridsize. 



3.4 The Lyapunov exponent for the sine map 

Loosely speaking, the Lyapunov exponent for a flow is the flow's mean rate of strain. Here we 
deflne the Lyapuov exponent exactly and compute it for the special case of the two-dimensional 
sine flow. We provide an asymptotic expressions for the exponent, in the limit of small and large 
stirring amplitudes. 



The integral of the period-one sine flow (3.14) gives a map from the 2-torus to itself: 



Xn+i = Xn + aa sin + 0„) , 

(3.15) 

Un+l =yn + Ot Sm [Xn+l + V'n) • 



By studying this map we may compute the Lyapunov exponent of the flow (3.14). Here a and 
a are stirring parameters and 0„ and ipn are phases that are randomized once per period. The 
chaoticity of this map is quantified by the function 7^: 

7, (a, a) = hm - log ( ), g G [0, 00) , (3.16) 

where w = {xq, i/q) is a two- vector of initial conditions, a is the stirring parameter and J™™ is the 
matrix 

j-cum/ ^ ^ {dxn+i,dyn+i) {dx^, Ovn) {dxi,dyi) 
{dxn,dyn) {dxn~i,dyn-i)"' {dxo,dyo)' 

with (dxiydyi) / {dxo,dyo) = l2x2- In what follows, for brevity, let 

{dXn,Oyn) W^W 

Finally, let the angle brackets (..) denote an ensemble average over all possible phase angles. 
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We study the function 7^ for particular values of q. For example, by setting g = we arrive at an 
expression of type 0/0, and applying L'Hopital's rule we replace - log(L''/^) with 

Thus, 

7o(a,a) = A= lim — (logL„); 72 (a, «) = lim — log(L„). (3.19) 

n— >oo Zn n— >oo ZTl 

The quantity A is the Lyapunov exponent of the map. In general, 7^ is a monotone increasing 
function of q (see for example, Finn and Ott [91]), so that 72 is an upper bound for A. This 
nonequality also indicates the importance of the order in which the averaging is performed. 




10 ' 10" 10' 10 ' 10" 10' 0.5 J 1.5 2 

loga ]oga a 

(a) (b) (c) 

Figure 3.1: (a) Lyapunov exponent 70 and its small- and large-amplitude forms; (b) Comparison 
between the Lyapunov exponent and 72; (c) Comparison between the Lyapunov exponent for the 
sine map with and without the randomization of phases. 



It is straightforward to compute 72: we obtain the exact result 



72 = 2 log 



l + ^ + ^yi6+^ 



(3.20) 



where f3 = a^fa. The evaluation of 70 is not so straighforward and so we have recourse to 
asymptotics. Thus, we replace the matrix J„ by its large-amplitude approximation 











ao?AnBn 



where = cos{yn + (pn) and = cos(a;n,+i + ipn)- The quantity L„ is easily calculated: it is 



Vo 



The order in which the computation is carried out is important: we take the log first and then 
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ensemble average over all phase angles. Thus, 

1 1 " 1 ^ 

- log L„ = log «^ + log a + - g log AlBl + - log ^ . 

It remains to compute the phase average (Y^^^^ log A^B^) = {log AlBf) + ... + {log A'^B'^). We 
start with the last term and integrate over 0„ and ^pn first, since this is the only term depending 
on these variables. Hence, 

{log AIBI) = — / d(f)n / dlpn log [COS^ (?/„ + 0„) COS^ {Xn+1 + Ipn)] , 
471" Jo Jo 

or 



(logA^i?^,) 

1 

~ 4^ 



27r /'27r p27r 

27r / logCOS^ {l/n + 0n) + / #n / dlpn log COS^ {Xn+1 (0n) + ^n) 

Jo Jo 



We use the integraQ 



27r 

2 



(iw log COS u = —An log 2 







to obtain the average (logA^i?^) = —4 log 2, so that J2^=i ^ogA'^Bf) = —2 log 2, and hence the 
result 

A~loga^ + loga-log4, a>l (3.21) 
The calculation in the limit a — > is not as simple, and so we resort to numerical methods. 



We compute 72 numerically and compare the result with the theoretical prediction Eq. (3.20). 
Hereafter a = 1. We find that 72 fiuctuates noticeably compared with 70. It is systematically below 
the theoretical value. However, by increasing the number of realizations in the ensemble average. 



this difference disappears. These results are presented in Fig. |3.1[ The reasonable agreement of 
the theory curve with the numerical one confirms the usefulness of the numerical calculation. Of 
significance is the fact that 70 is nonnegative even for very small a. Thus, the sine map is chaotic 
for all parameter values. 

The large-a result for 70 is checked against a numerical simulation and a perfect agreement is 
obtained. In Section 13.51 we shall need a formula for the strain rate for moderate values of the 
stirring amplitude a and so we fit a curve to 70 (a) for the parameter range a G [0, 1]. To make the 
system (1) invariant under a —>■ —a, the fitted curve should be a polynomial in a^. The equation 

A(a,a) = O.llSaa^ (3.22) 

is found to be optimal. Going to higher order in does not improve the accuracy of the fit. 

As a final computation we contrast the sine map with random phases with one in which the phases 



*For a derivation of this integral relation, see Appendix A. 
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are constant and set to (p = n, ip = 0. This contrast is shown in Fig. |3.1[ c). For large a, these maps 
have the same behaviour, growing as loga^, and this is explained by the analysis of this section. 
For small a there is an approximate parameter range a G [0, 0.5] for which the constant-phase map 
does not exhibit chaos. The asymptotic form A = lim^^oo (1/2?T') (logL„) is arrived at after 40 
iterations for the random-phase map and after 10,000 iterations for the constant-phase map. Thus, 
in contrast to the constant-phase map, the random-phase map is efficient in separating trajectories 
at all stirring amplitudes and over all but very short timescales. 



3.5 Numerical simulations 



We integrate Eq. (3.1) in the manner explained in Sec. 3.3 The initial conditions are chosen to 



represent the sudden cooling of the binary fluid: above a certain temperature Tc the homogeneous 
state c = is stable, while below this temperature the mixture free energy changes character to 
become | (c^ — 1)^, which makes the homogeneous state unstable. Thus, at temperatures T > Tc, 
c{x) = + [fluctuations], and the sudden cooling of the system below leaves the system in 
this state. In order for the CH equation without thermal noise (3.1) to describe the evolution 
of these initial conditions, the cooling we have imposed must take the system to T = 0. Thus, 
we are working in the so-called deep-quench limit. With the initial conditions c{x,0) = + 
[fluctuations], both components of the fluid are present in equal amounts and the spatial average 
of the concentration fleld is zero. 

Although the case without flow has been investigated before ^U\, we reproduce it for two 
reasons. First, it serves as a validation of our algorithm. Second, we shall conflrm the free 



energy laws (3.9) and (3.11) which to our knowledge have not been examined before. Finding an 



appropriate measure of bubble growth will be important for understanding the behaviour of the 
bubbles when we introduce stirring. 



The results of these simulations are presented in Fig. 3.2 



The scaling function kfs {k/ki,t) is 
approximately time independent for t > 20,000 as evidenced by Fig. 3.2[ b), implying that the 
scaling exponents for ki and F are to be extracted from the late-time data with t > 20, 000. For 
t > 30, 000 flnite-size effects spoil the scaling laws. Thus, while the scaling exponents are extracted 
from a small window, the flt is good and this suggests that the power laws obtained give the true 
time-dependence of ki and F at late times. In this way we recover the dynamical scaling regime 
in which kfs {k/ki,t) is time independent, in addition to the decay law ki ~ t~^^^ [5]. Running 
the simulation repeatedly for an ensemble of random initial conditions improves the accuracy of 
the exponent. Furthermore, we obtain the energy laws F ~ cit-^^^, and 1 - ~ cst"^/^. Thus, 



a 
~F 



C2 



^-1/3 



.1/3 



Cl 



for t > 1, 



(3.23) 



so that cT^/F, 1/F and 1/ki all grow as t^^^ and hence provide identical measures of scale growth, 
in agreement with the classical assumption that there is a unique lengthscale in the problem [U [5] . 



We investigate the stirred case in a similar manner and obtain the results below, after t = 30, 000 
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Figure 3.2: (a) Concentration, at t = 30,000; (b) Structure function kfs {k/ki,t) in the scaling 
regime t > 20, 000; (c) Time dependence of ki, with power law exponent close to the LS exponent 
1/3; (d) Free energy, F, exhibiting the time dependence F ~ t"*^'^^^, close to the LS exponent. 



timesteps. In order to ensure that we are in a steady state, we study the inverse lengths F and ki 
and the tracer variance (see Eq. (3.10)) which measures the homogeneity of the concentration 

^2 



1^0"' = for a homogeneous mixture) |24] |27] . These quantities are time dependent and have the 
property that after some transience, they fluctuate around a mean value, without secular trend. 
This can be seen in Fig. 3^ In Fig. 3.4[ for a = 0.1, the concentration field at late time looks 
similar to that for a = 0. This is because for small a, the effect of advection is diffusive, with small 
movements of particles giving rise to bubbles with jagged boundaries, but no breakup. The PDFs 
of the concentration for a G [0, 0.01, 0.2] have the same structure. Nevertheless, it is clear from 

[0.01, 



Fig. |3.3[ b) that coarsening arrest takes place for a = [0.01, 0.2], in contrast to a = 0. As a 
increases, the bubbles become less and less evident. For large a, we see regions containing filaments 
of similar concentration. After a period of transience, the free energy F, the mean wavenumber 
ki and the variance cx^ fluctuate around mean values, confirming the existence of the steady state. 
The bubble size is extracted from the quantity (a'^/F), where (■) denotes a time average over the 
fluctuations. We investigate how the mean lengths scale with the Lyapunov exponent A. 



In Fig. 3.5 we see that for small A the proxy (a'^/F) for bubble radius scales with the LS exponent 
as A"^/^, suggesting an equilibrium bubble size on this scale. For larger A, the hyperdiffusion 
breaks up these bubbles and mixes the fluid. This process leads to a faster-than-algebraic decay of 
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bubble size, visible in the figure. Other measures of the bubble size, such as the mean wavenumber 
{ki) or the quantity 1 — (cr^) do not possess a clear scaling law. For example, in Fig. 3.5[ b) we 
plot log(/ci) against log A. The exponent changes over the three decades of data, indicating that 
the wavenumber (ki) does not have a clear scaling law with the Lyapunov exponent A. However, 
the dependence of the time-averaged free energy on the Lyapunov exponent has the power law 
behaviour (F) ~ A°'^^ over three decades, suggesting a genuine power law relationship. The time 

average of the Batchelor wavenumber kf, = [| Vc|^/c^]^/^ [52] also possesses a clean power law over 
three decades. 



In Sec. 2.2 we showed that in a situation of dynamical equilibrium characterized by a single 



lengthscale, the potential part of the free energy Fp = | J d'^x (c^ — 1) should be in equipartition 
with the kinetic part = ^ J d?x^ |Vc|^. In Fig. 



3.5 



we see that this is not the case when the 
his is an indication of a crossover between 



bubbles start to break up as the stirring is increased 
the bubbly regime and the well-mixed one in which the variance is reduced by stirring 



In Fig. 



3.5 



the breakdown in the power law relationship (a'^/F) ~ X^^^^ for large stirring amplitudes 
suggests a crossover between bubbly and diffusive regimes. Another way of seeing the transition 
between these regimes is to study the stationary probability distribution function (PDF) of the field 



c {x, t). We do this in Fig. 3.6, For small values of a we note that the PDF has sharp peaks at ±1, 
indicating the effectiveness of phase separation at these stirring amplitudes. For a = 1 however, 
the PDF has a Gaussian core, indicating that genuine mixing by advection-hyperdiffusion is taking 
place. For intermediate values of a the PDF is a combination of these two different distributions 
of concentration. 



3.6 A one-dimensional equilibrium problem 

Having observed the crossover between the bubbly and the well-mixed regimes, we study a one- 
dimensional model to shed further light on the process by which this occurs. 



We examine the archetypal hyperbolic flow v = {—Xx,Xy), with strain rate A. This flow tends 
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(d) (e) 

Figure 3.4: The steady-state concentration field after t = 30, 000 timesteps, for various stirring 
amplitudes a. 

to homogenize the concentration in the y-direction through stretching. At late times the problem 
then becomes one dimensional: 



Scaling lengths by ^7 and restricting to the steady case, Eq. (3.24) becomes 



/ 7 \ ac a / q N a c 
- A a;— = — (c^ -c) - — . 3.25 



There are two sets of boundary conditions (BCs) that are of interest. 



Bubble BCs: c (±00) = -1 + 5, (±00) = 0, and 5 > 0. 
Front BCs: c(±oo) = ±1, and c' (±00) = 0. 



The front solution has the interpretation of being the transition layer between filaments, the width 
of the layer giving the width of the filament. Before carrying out numerical work, we shall examine 
two asymptotic regimes where analytical progress is possible: the case without stirring, and that 
without phase separation. 
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Figure 3.5: (a) The plot of log((T^/F) against log A suggests the scaling law (a'^/F) ~ A^^/'^ for small 
stirring amplitudes while for large amplitudes this proxy for the bubble radius decays exponentially 
with increasing Lyapunov exponent, (b) Graph of log(/ci) against log A, with approximate power 
law behaviour (ki) ~ A°'^^. The exponent is not the same over the entire range, however, (c) 
The graph of log(F) exhibits a much cleaner power law behaviour (F) ~ A°'^^. (d) The Batchelor 

wavenumber ki, = [|Vc|^/c^^^^ also possesses a clear power law behaviour {ki,) ~ A°-^^. 



No stirring 



In the absence of advection (A = 0) the model equation reduces to 



which integrates to 



c — c — 



dx'^ 



0, 



Ax + 5, 



(3.26) 



(3.27) 



where A and B are constants to be determined by the boundary conditions. We set bubble BCs 
c(±c>o) = — 1 + 5, and (±oo) = 0. Using the bubble BCs and the requirement that the second 
derivative of c should not diverge at infinity, A is set to zero. Thus we have the Newton equation 



dx"^ 



d_ 
dc 



-l{c'-lY + Bc 



(3.28) 
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a=0.1 a=0.3 a=0.5 




c c Id 



(d) (e) (f) 

Figure 3.6: Normalized PDF of concentration in the steady state, (a), (b) Segregation-dominated 
flow (a = 0.1,0.3); (c), (d) Crossover to quasi- diffusive regime (a = 0.5,0.7); (e) Quasi-diffusive 
regime {a = 1.0); (f) Semilog plot of PDF for a = 1.0. Note Gaussian core and exponential decay 
down to the cutoff at |c| < 1. 



By multiplying both sides by dc/dx we obtain the energy relation 




(3.29) 



where we identify the potential energy V (c) = — ^ (c^ — 1) + Be and the total energy E. The 
bubble boundary conditions then fix the values of E and B. The bubble or homoclinic solution is 
obtained when two roots of the equation E = V (c) coincide with one extremum of V, given by 

-c = B. 

We pick out the homoclinic orbit by letting this extremum be at c = — 1 + 5, giving B^ = 
(— 1 + 5)^ — (1 — 5). In addition, we specify the energy to be 



E = V{-l + 6) 
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Then in an interval (ci = (— 1 + 5) , C2), indicated in Fig. 3.7 (a), the energy E exceeds V (c), and 
E — V > 0, with equahty at the interval boundary. Thus, the homoclinic solution is given by 



X — Xj- = ± 



V2 



dc 



V{-l + 5)-Bsc+\{c^-l) 



By taking the reference level Xj. = 0, we obtain the equation 



x = ± 



V2 



c{x) 



dc 



(0) JV{-l + 5)-B,c+\{c'-iy 



(3.30) 



Equation 3.30 specifies c(x), which can then be obtained numerically, as in Fig. 3.7 (b). In the 




Figure 3.7: (a) The potential function V {c) = — | (c^ — 1)^ + Be. Here we set B = (— 1 + 5)^ — 



(— 1 + 5). Notice the well extremities at ci = — 1 + 5 and C2; (b) Integration of Eq. (3.30) for 
5 = 0.001 giving the bubble profile (homoclinic orbit); (c) Periodic solution for E < V{ci), 
Ci < c < C2, which is close to the homoclinic orbit, (d) Phase space diagram of solutions: the 
homoclinic and large-amplitude periodic solutions nearly coincide. A small-amplitude periodic 
solution is included for reference. In the homoclinic case, the zeros of Cx lie at the turning points 
ci and C2. 

language of dynamical systems, the solution is a homoclinic orbit: the system returns to its initial 
state at X = —00 after infinite time. The phase portrait of the solution is presented in Fig. |3.7Kd). 
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Finally we note that Eq. (3.28) possesses a host of different solutions, depending on the boundary 
conditions imposed. Let us take B as before but leave the energy unspecified. We study the 
Newton equation c" (x) = —dV/dc, where V (c) is the potential V (c) = — | (c ^ — 1)^ + Bsc. 
Periodic solutions exist when the system is placed in the potential well of Fig. |3.7| (a), and when 
E < V {—1 + 6) . To obtain a solution with the fiat peaks characteristic of a bubble, the system 
must access the anharmonic part of the potential, so we take E < V {—1 + 6). Thus the phase 
portrait of this periodic bubble solution almost coincides with the homoclinic orbit. On the other 
hand, for energies close to the well minimum V (cq) = Vq the periodic solution is a small sinusoidal 
undulation. These results are best summarized by introducing the period of oscillation, given 
by iH 



where ci and C2 are those roots of the equation V (c) — E = indicated in Fig. 3.7 



For E >Vq the oscillations are small and the period is l-n^V" (cq). 

For E <V {—W S) the solution is still periodic with energy-dependent period. This is the 
one-dimensional gas of Fig. 3.7[ c). 



For 1 + 5), T— > +oo, giving the homoclinic solution. 

For 5 — i> and E (^—\ ^ 6\ T— > —ioo, which corresponds to the heteroclinic or front 
solution. 



A quasi two-dimensional problem 

Following the mechanical intuition developed in the previous section, we outline briefly what 
happens in the quasi one-dimensional setting when the full two-dimensional problem has azimuthal 
symmetry. In this case, the stationary CH equation is 

[c'^-c- Vjc] = 0, 

where is the radial part of the Laplacian, 

^ r dr Or ^ ^ 

Using the same analysis as before, we obtain the equation — c — V^c = B, which is a Newton 
equation with damping: 

d?c ^Idc dV 
dr"^ r dr 9c ' 

where V (c) = — | (c^ — 1)^ + Be is the potential. The large-amplitude and homoclinic solutions 
presented previously are modified by the presence of damping: at small r there is significant 
dissipation proportional to 1/r while the far-field problem is one of small oscillations with damping. 
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Figure 3.8: (a) Phase portrait for the formerly homochnic case. The solution spirals towards the 
attractor (cq, 0); (b) Profile of the solution and its derivative; (c) The full two-dimensional solution 
with azimuthal symmetry. 

These oscillations have period I'k^/V" (cq) (cq is the well minimum), and decay to the attractor 
(co, 0). These results are shown succinctly in the phase portrait of Fig. 3.8[ a). 



The no-segregation limit 

By measuring lengths in terms of the diffusion lengthscale (7D/A)^^^, Eq. (3.25) becomes 



which in the limit of large A reduces to 



XV 



(3.33) 



(3.34) 
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where c (x) = c (0) + dx'v (x'). Equation (3.34) is a generalized Airy equation with solution [93] 



v{x) = ^ CyVy {x) , (3.35) 
where the functions Vy (x) are given by 

/■oo 

Vu{x) = ei, I exp (e,^xt — 1^/4) (it; z/ = 0,1,2,3. (3.36) 
Jq 

The phases are the fourth roots of unity: = exp {TTh'i/2) and the constants Cu must satisfy 
^j^C^ = 0. A plot of vi and ^2 is shown in Fig. 3.9 By studying the phases e,y, we see that vi 
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Figure 3.9: (a) Generalized Airy function Vi (x), with real and imaginary parts shown; (b) The 
function V2 (x), which is pure real. 

and V2 are related to and vq respectively by a parity transformation. 

We set Ci = —C3 = 1 and Cq = C2 = to obtain the front solution of the hyperdiffusion equation: 



Cfront (•^) 



[vi (x') — V3 (x')] dx' . 



(3.37) 







Simplifying and restoring the dimensional units, this is 

r>00 



Cfront (x) 



t 



sin (tx) exp 



4A/7D 



dt. 



(3.38) 



By inspection of the phases eu, it is clear that the hyperdiffusion equation has no bubble solutions 
on the whole real line. In Fig. 3.10| we present front solution. 



Finally we note that the generalized Airy functions Vi, (x) give the solution of the nonstationary 
advection-hyperdiffusion equation 

dc dc d^c 



dt ^ dx 



dx^ 



where we have again expressed lengths in terms of the diffusion scale {^D /Xf'^^ . Time is expressed 
in terms of the Lyapunov time 1/A. By linearity, we may separate the variables, writing c(x,t) = 
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Figure 3.10: (a) Front solution (3.38) obtained through quadrature 



e *c (x). The Airy functions (x) = e'^'^^'^-'rft satisfy 



dx^ dx 



dt 



= V, {x) , 



(3.39) 



which is the equation for c{x). Since vq and V2 do not have the correct far- field behaviour, we are 
restricted to the dimensional solutions 



c (x, t) 



X 



(3.40) 



The functions f 1^3 (x) have the asymptotic form t>i 3 (x) ~ — 1/x + 6z/x'', for |x| ^ 1 [51], and 
therefore have the correct far-field behaviour. Thus, in the absence of segregation, a bubble profile 
will decay in time owing to the twin processes of advection and hyperdiffusion. 



Numerical simulations of the boundary value problem 



Armed with this knowledge, we attempt a numerical solution of the boundary- value problem (BVP) 

- A ( ^ ) x^ = ^ fc' -c) - (3.41) 



7 \ dc 



dx 



dx^'' 



where lengths are measured in terms of the transition layer width The analytic study of 

the hyper diffusive limit suggested that seeking a bubble solution may not be sensible, and indeed 
it is not: solving the time-dependent version of Eq. (3.41) with bubble initial conditions, the 
concentration profile always decays to a constant. We therefore focus on the frontal boundary 
conditions. 



We solve the BVP (3.41) with frontal boundary conditions in Fig. 3.11 The front steepens with 
increasing X'j/D. Physically, this corresponds to a strengthening of local shear in the velocity, 
which forces the filaments of fluid to become narrower. The numerical solution agrees exactly 
with the analytical solutions in the limit of large and small A. At large A, the front width has 
the hyperdiffusive scaling, indicating a regime dominated by hyperdiffusion, in agreement with 
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the full two-dimensional numerical simulations in Sec. 3.5 We note that the concentration profile 




o 10" 




(a) 



(b) 



Figure 3.11: (a) Front profile for \-f/D = 0.001,0.01,0.1,1,10,100,1000. The front steepens as 
X'y/D increases; (b) Front "width as a function of the parameter X'-f/D: for large parameter values, 
the front width scales hyperdiffusively as {X'f /D)^^'^. 



takes values |c| > 1 in some places, a consequence of the fact that Eq. (3.25) has no maximum 



principle, in contrast (say) to the advection-diffusion equation. Again, this result is consistent with 



the numerical results of Section 3.5, where we find maxg-go |c(a;,t)| > 1. This superabundance of 



one binary fluid element is unreasonable on physical grounds, pointing to a shortcoming in the 
advective CH model. However, we shall see in the next section that the qualitative features of the 
constant-mobility model and those of the more realistic variable-mobility model are in agreement, 
suggesting that the simpler model gives an adequate description of advective phase-separation 
dynamics. 



3.7 The variable- mobility case 

To gain further insight into the dependence of the bubble size on the Lifshitz-Slyozov and Lyapunov 
exponents, we study the variable-mobility Cahn-Hilliard equation 

dc 

_ + z;.Vc = V-[D(c)V/i], (3.42) 

where as before /i = — c — 7V^c. By introducing this additional feature, we recover a system 
that, at least in the case without flow, has a solution c{x,t) confined to the range [—1, 1] [8]. 

We specify the functional form of the mobility as 

D (c) = Do{l- vc') , (3.43) 

fixing 1] = 1. This modification corresponds to interface-driven coarsening because inside bubbles 
(c = ±1) the mobility is zero. For other values of t] we have a mixture of bulk- and interface- 
driven coarsening. This equation has been investigated by Bray and Emmott |l96J for v = 0, in an 
evaporation-condensation picture. They show that the typical bubble size i?b (t) grows as t^^^, a 
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result limited to dimensions greater than two. Numerical simulations [5] suggest that this growth 
law also holds in two dimensions. We couple the modified segregation dynamics to the sine fiow. 

For i; = 0, we obtain the growth law for the scale Rh = l/^i, namely i?b ~ The growth 

exponent is close to the value 1/4 predicted by the LS theory. In addition, we obtain the free 
energy decay laws F ~ ^-o-267 ^gain close to the LS value. Thus, as in the constant-mobility case, 
a description of the problem lengthscales can be given in terms of the bubble energy or in terms 
of i?b, with identical results. 

Because the variable-mobility calculation is computationally slower, we shall perform the numerical 
experiments with fiow at lower resolution. We do not anticipate that this will change the results, 
because in integrations of the constant-mobility equations at resolution 256^ the scaling exponents 
were unaffected. Switching on the chaotic flow, we observe the steady state, characterized by 
the fluctuation of the free energy, the wavenumber ki, and the variance around mean values. 




Figure 3.12: (a) The plot of log(cr^/F) against log A suggests the scaling law (a'^/F) ~ A"^/"^ 
for small stirring amplitudes, while for large amplitudes this proxy for the bubble radius decays 
exponentially with increasing Lyapunov exponent, (b) The graph of log(F) exhibits power law 
behaviour (F) ~ A'^'^^. The predominance of kinetic free energy over potential free energy at large 
A is visible. 



As in the constant-mobility case, the graph of the time-averaged free energy and time-averaged 
Batchelor scale possess clear scaling laws, while the graph of the mean wavenumber (ki) does not. 
For small A, the proxy for bubble radius (a'^/F) scales approximately as A~^/^, suggesting a decay 
of bubble radius according to the LS exponent. Thus, the analogous result of Section [3l5] is not 



fortuitous. For large A this quantity decays exponentially to zero at a rate (with respect to A) 
close to that of the constant-mobility system. In each case, the flt of the data is not as good as 
the constant-mobility case, the computational error being larger here because we have run the 



simulations at resolution 256^. We present these results in Fig. 3.12 



The results are an exact replica of those for the constant-mobility case. The arrest of bubble growth 
for small A is therefore a genuine phenomenon whose behaviour depends on the LS exponent of the 
corresponding unstirred dynamics. The mixing at large A is identical to that for which the mobility 
is constant. Thus, while the variable-mobility model is more realistic than its constant-mobility 
counterpart, their properties with regard to chaotic stirring and eventual mixing are the same. 



3.8. Summary 



68 



3.8 Summary 

We have shown how, in the presence of an external chaotic stirring, the coarsening dynamics of 
the Cahn-Hilhard equation are arrested and bubbles form on a particular lengthscale. The system 
reaches a steady state, characterized by the fluctuation of the free energy F, the variance a^, and 
the wavenumber ki around mean values. A measure of the typical bubble size is given by the time 
average of a'^/F. For sufficiently large stirring intensities, the bubbles diminish in prominence, 
to be replaced by filament structures. These are due to the combination of the advection and 
hyperdiffusion terms in the CH demixing mechanism. In this regime, the concentration tends to 
homogenize, so that the stirring stabilizes the previously unstable mixed state. 

We have explained with a one- dimensional model how this transition arises and how the filament 
width in the homogeneous regime depends on the Lyapunov exponent of the chaotic fiow. This 
latter investigation exhibits the local superabundance of one of the binary fiuid components, which 
we reject as unreasonable. In the case of no fiow this difficulty was overcome by the addition of 
a variable mobility [S]. However, the quasi-diffusive regime identified in our simulations is robust 
in the sense that it is present in both the constant- and the more realistic variable-mobility cases. 
We therefore expect to see this remixing phenomenon in real binary fluids. 



Chapter 4 



Estimating mixedness in a stirred 
Cahn— Hilliard fluid 

4.1 Overview 

In this chapter we study the advective Cahn-Hilhard equation in the passive setting, in the presence 
of sources and sinks of matter. We focus on the symmetric mixture, which contains equal amounts 
of each binary fluid component. We have mentioned the often undesirable coarsening tendency 
of the binary fluid, and the efforts made to suppress it. We therefore introduce a quantitative 
measure of coarsening suppression by studying the p^^ power-mean fluctuation of the concentration 
about its average value. By fluctuations about the average value, we mean spatial fluctuations 
around the average spatial concentration (which is zero for a symmetric mixture); we then average 
these fluctuations over space and time. In effect, we study the time-averaged norm of the 
concentration. If this quantity is small, the average deviation of the system about the well-mixed 
state c = is small, and we therefore use this quantity as a proxy for the level of mixedness of the 
fluid. 

In many applications, the Cahn-Hilliard equation is driven by some external effect, for example, 
chemically patterned substrates ^7\, electrically charged media [HS], or the subject of this report, 
stirring. In this chapter, we focus on phase separation in the presence of advection and a source 
term in the Cahn-Hilliard equation. This source term arises in a number of ways. In [88] it is 
maintained by thermal diffusion, through the Ludwig-Soret effect [99], in which concentration 
gradients are induced by imposed temperature gradients. One can also produce concentration 
gradients by simply injecting matter into the system. The approach we take here is motivated 
by studies of the advection-diffusion equation [281 EHl |30l EI]. There, the problem is a linear one, 
and describes miscible liquids, while fluctuations about the mean are measured by the variance or 
centred second power- mean of the fluid concentration [32l [33] . The variance is reduced by certain 
stirring mechanisms. By specifying a source term, it is possible to state the maximum amount by 
which a given flow can reduce the variance, and hence mix the fluid. By quantifying the variance 
reduction, we can classify flows according to how effective they are at mixing. 
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4.2 Measures of mixedness 



In this section, we introduce the advective Cahn-Hilhard equation with sources and recall some the 
properties of this equation. We outline the tools and notation we shall use to analyze concentration 
fluctuations. The results will be valid in any spatial dimension. 

We study the advective Cahn-Hilliard equation in the presence of flow, and for prescribed sources 
and sinks, 

dc 

— + v\/c = DA{c^-c--fAc)+s{x). (4.1) 

Here D is Cahn-Hilliard diffusion and is the typical thickness of transition zones between 
phase-separated regions of the binary fluid. 



In this chapter we work with a nondimensionalization of Eq. (4.1) that leaves three control pa- 
rameters in the problem. Therefore, we can unambiguously study limits where control parameters 
take large or small values. Let Tq be a timescale associated with the velocity v {x,t), and let Vq 
be the magnitude of v {x,t). Let 5*0 be the magnitude of the source variations and finally, let L 
be a lengthscale in the problem; for example, if the problem is solved in a hypercube with peri- 
odic boundary conditions, we take the lengthscale L to be the cube length. It is then possible to 



write down Eq. (4.1) using a nondimensional time t' = t/T^ and a nondimensional spatial variable 
x' = x/L, 

^ + v^v ■ Vc = D'A' [c'-c- VA'c) + S',~s {x') , 

V and s are dimensionless shape functions, Vq = TqVq/L, D' = DTq/L"^, 7' = 7/L^, and where 
5*0 = SqTq. As before, the quantity D' = DTq/L'^ is identified with the ratio Tq/Td, being the 
ratio of the velocity timescale to the diffusion timescale. Following standard practice, we shall now 
work with the dimensionless version of the equation, and omit the prime notation. For ease of 
notation, we shall henceforth take v to mean V^v. 



We recall the following properties of (4.1) from Ch. [2j These will be useful in what follows. 
• If the source s [x) is chosen to have spatial mean zero, then the total mass is conserved, 

[ c{x,t)rx= I DA{c^ -c-jAc)rx+ I s{x)rx 
Jn Jn Jn 



d 
dt 



+ [boundary terms] 



where VL is the problem domain in n dimensions and \VL\ is its volume. The boundary terms 
in this equation vanish on choosing no-fiux boundary conditions n ■ Vc = n ■ V/i = on dVL, 
or periodic boundary conditions. Here n ■ V is the outward normal derivative. For simplicity 
we consider the periodic case, although this is not necessary. 

There is a free-energy functional 

F[c]= f \\{c^ -if + \^\Wcf\d''x, /i= ^ =c=^-c-7Ac, (4.2) 

In L J OC 
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where fi is the chemical potential of the system. For a smooth concentration field c{x,t), 
the free energy satisfies the evolution equation 

^ dF 



dt 



-D / |V/xr(i"x+ / /i (-V ■ Vc + s) c/"x, 
'n Jn 



and in the absence of sources and stirring decays in time. 

To study the spatial fluctuations in concentration, we consider the power means of the quantity 
c {x, t) — \^\~^ Jq c {x, t)d^x, 



Mp{t) 



c(x,t) — — / c{x,t)d"'x 
l"l Jn 



(4.3) 



For a symmetric mixture in which c {x, t) d^x = 0, this is simply 



\c{x,t)fd''x 



where we have introduced the norm of the concentration, ||c||p. The quantity Mp is a measure 
of the magnitude of spatial fluctuations in the concentration about the mean, at a given time. 
Since we are interested in the ultimate state of the system, we study the long-time average of 
concentration fluctuations. We therefore focus on the power-mean fluctuations 



where (■) is the long-time average 

{■) = }imj[\.)ds, 

t^OO t Jq 

provided the limit exists. We shall repeatedly use the following results for the monotonicity of 
norms, 



1 1 

< \n\p i 



9' 



\gis)\ds < 



\g {s)\''ds 



i<P<g,/eL« (fi), 

q>l,geL'^i[0,t]) 



(4.4) 



discussed in Ch. [2l 



The Cahn-Hilliard equation and its free energy functional contain high powers of the concentration 
c (c^ and respectively), and we can therefore estimate rrip for speciflc values. In particular, in 
the following sections, we shall prove the following result in n dimensions: 



Given a smooth solution to the advective Cahn-Hilliard equation, the long-time average 
of the free energy exists, and therefore rrip exists, forp G [1,4]. 
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The only constraints we impose on the external forces are that the velocity, its derivatives, and the 
source term be always and everywhere bounded. We take our result one step further by explicitly 
evaluating upper and lower bounds for m^, and this gives a way of quantifying concentration 
fluctuations in the stirred binary fluid. 



4.3 Existence of long-time averages 

In this section, we prove a result concerning the existence of the long-time average of the free 
energy, and of the power means rUp, for p G [1, 4]. 



Given the velocity field v{x,t) G H '°° (il), the source s (x) G L {^l), and smooth 
initial data for the advective Cahn-Hilliard (4.1), the long-time average of the free 
energy exists, and thus rrip exists, forp G [1,4]. 



The proof relies on the free-energy evolution equation. Using this law, we find uniform bounds on 
the finite-time means {F)t and (M|)i, where 

{■)t=j [\-)ds, (■)= lim(-),. 

Using the monotonicity of norms, the uniform boundedness of {M^)t, follows, for p G [1,4]. The 
proof proceeds in multiple steps, which we outline below. 



Step 1: Analysis of the free-energy evolution equation 

Given the conditions on the external forces outlined above, and the results in Ch. [2| there is a 
unique smooth solution c{x,t) to the advective Cahn-Hilliard equation (4.1), at least for finite 



times. Thus, we turn to the question of the long-time behaviour of solutions. We exploit the 
smoothness property of the concentration field c{x,t) in formulating an evolution equation for the 
free energy. The free energy of the Cahn-Hilliard fiuid is 



F[c] 



{^-lf + \^\Vc\ 



1 /'^2 
4 



d'^X. 



Given the smooth, finite-time solution c(a;,t), we differentiate the functional F [c] with respect to 
time and obtain the relation 



dF 
~dt 



-D [ |V/i|^c/"x+ [ /i (-V ■ Vc + s) 
Jq Jn 
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using the no-flux or periodic boundary conditions. By averaging this equation over finite times, 
we obtain the identity 



{F)t + D{ / |V/irrf"x 



fiscT'x 
n I t 



n I t 



(4.5) 



We single out the quantity (F) for study. Owing to the nonnegativity of F{t), we have the 
inequality (F) > 0. Therefore, we need only consider two possible cases: (F) = 0, and (F) > 0. 
We shall show that (F) 7^ is not possible, and in doing so, we shall produce a uniform {t- 
independent) upper bound on {F)t. 

Let us assume for contradiction that (F) > 0. Then, given any e in the range < e < (F), there is 
a time Tg such that (F) —£< {F)t < (F) +e, for all times t > T^. Thus, for times t > T^, the time 
average {F)t is strictly positive. Henceforth, the inequality t > is assumed. Using V ■ i> = 0, 



Eq. (4.5) becomes 



{F)t + D{ / |V/i|'rf"x 



(4.6) 



As in Sec. 2.6, we use the following string of relations. 



Acv ■ Vcd^'x 



(didic) {vjdjc) d^x, 



(dic) [di (vjdjc)] d^x, 



{d,c){d,v,){d,c)d^x- / {d,c){vV){d,c)d''x, 
'n Jn 

= - wWw^d'^x, w = Vc, Wij = \ {diVj + djVi) 
Jq 

where we have used the antisymmetry of the operator v • V to deduce that 



bvV ■ M^'x = 0, 



for any function {x,t). The quadratic form wWw^ satisfies 

< nmax |iyjj|||iu||2. 



which gives rise to the inequality 



Acv ■ Vcd'^x 



< n \ sup \ Wij 



Wcrd^'x. 



(4.7) 



The matrix W is the stretching tensor |lUUj . the appearance of which highlights the importance 
of shear and stretching in the development of the morphology of the concentration field [20] . 
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For each time t' G [0, t], we split the chemical potential fi into a part with mean zero, and a mean 
component: ^ = Jl {t') + fi' {x, t'), where /i' {x, t') dJ^x = 0. Then, for any function {x, t) with 
spatial mean zero, we have the relation (f)fid"'x = (p^'dP'x. Using this projection relation, 
Eq. (4.6) becomes 



{F)t + D^j \Vii't(rx 
Owing to the positivity of (-F)t, we have the inequality 



li'scTx ) +{-f Acv- VccTx 

t \ Jo It 



D{ I |V/ifrf"x) <( / li'sdTx) +{l /^cv-VcdTx). 
In / 1 \Jn / 1 \ Jn It 



Finally, we employ the Poincare inequality for mean-zero functions on a periodic domain Vl 
[0,L]^: 



Combining Eqs. (4.7), (4.8), and (4.9) gives the following inequality: 



(4.9) 



(4.10) 



where W^o = ^'^Y>t&[o,oa),n,i,j There are no angle brackets around the source term because 

s [x) is independent of time. 



Step 2: Obtaining a bound on (||/u'||2)t 

Using the string of relations 

f 7 I Vc|^ (Tx = [ [fic + c^- c^] (Tx = I [/x'c + - c^] (Tx < 
Jn Jn Jn 

we obtain the inequahty 

7|Vc|^rf"x < ||^'||2||c||4+ \\c\\l. 



2 C 2 + C 



Combining Eqs. (4.10) and (4.11) gives rise to the inequality 



D 



27V 

T 



f^'\\'i)t<{\\f^'\\l)! \\s\\2 + n\n\iWMt)t +2\Qr^w^{\\c\\t) 



1 

4\ 4 



25 



(4.11) 



|4\ 2 
It 5 



4.3. Existence of long-time averages 



75 



1 

a quadratic inequality in (||/i'||2)t' • Hence, 



1 f L\' 



4\4 

t 



1 f L\ 



IN 2 



+ ^(^) W(||s||2 + n|i]|^w^oo(||c||^)|j>8D|fi|M^j w^oodlcllDl. 



'27r\ 



A less sharp bound is given by 

n\n\'^w^{\\c\\t)!) + ^ . ^ , 

which is an upper bound for (H/x'UDt, in terms of the forcing parameters, and (||c|||) 



1x2 8|^]|^ fL\ 



WMl)!, (4.12) 



Step 3: An upper bound on 



We have the free energy 



F\c\ 



{c'-lY + l^\Vc\ 



[W-\c'] rx + \\n\ 



[ [lcfi'{x,t')-lc'] rx + l\Q\>0. 
Jn 



Hence, 



/ c^rx <2 [ cfi'rx + \n\ < 2||c||2||/x'||2 + M . 

Jn Jn 



Time averaging both sides and using the monotonicity of norms (4.4), we obtain the result 



\c\\t)t<\n\ + 2\n\H\\c\\t)t{M)!- 



Using the bound for 

{\\c\\i)t < M 



/||2\ 



)t in (4.12), this becomes 



D \2tt J 



\s\\2 + n\n\^^Woo{\\c\\i)f ) +4nD\n\^2 



1 /27r 



L 



W^,(\\c 



|4\ 2 



We therefore have a t-independent equation for the upper bound on (||c||4)( 



\c\\t)t<mr^[v,s,D] 



(4.13) 
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where m™^^ solves the polynomial 



1 2 1 Z' 2 \ 

|s||2 + n|0|3 ^00^4°''''') +4nL»|17|5 f -^J TVoo(m4 



2 

(4.14) 



The highest power of m^^^ on the left-hand side is (m™^^)^, while the highest power of m™^^ on 

3 

the right-hand side is (m™'^^) ^ . Thus, this equation always has a positive solution, and moreover 
this positive solution is unique. 

We obtain the following chain of uniform (t-independent) bounds. Each bound follows from the 



previous bounds in the chain, and the first bound follows from Eq. (4.13) 



(||c||4)t is uniformly bounded, 
(llcllDt is uniformly bounded, 
(||/i'||2)t is uniformly bounded, 
(||Vc||2)i is uniformly bounded, 
{F)t is uniformly bounded. 



(4.15) 



for all t > Tfr. Owing to the uniformity of these bounds, they hold in the limit t —>■ oo. The 
result (F) < oo implies the existence of a uniform bound for F{t), almost everywhere. Given 
the differentiability of F{t), this implies that F (t) is everywhere uniformly bounded, and thus, 
(F) = 0, which is a contradiction. Therefore, the only possibility for (F) is that it be zero. It is 
straightforward to verify that by taking {F) = 0, and making slight alterations in Steps 1-3, the 



bounds in Eq. (4.15) still hold. 



Let us examine the significance of our result. We have shown that for sufficiently regular flows 
and source terms (specifically, v {x,t) G H^'°° (Q), for all t G [0,cx3), and s (x) G (fi)), there 
is an a priori bound on the free energy {F [c]). We have shown that the system always reaches a 
steady state, in the sense that {F) = 0. We have also found an upper bound for the measure of 



concentration fluctuations, as the unique positive root of the polynomial equation Eq. (4.14). This 
bound depends only on the source amplitude, the diffusion constant, and the maximum stretching 
Woo. Using the monotonicity of norms, this number serves also as an upper bound on rup for 
p G [1,4]. Let us comment briefly on the volume term in the equation (m™'^^)*^ = + .... Since 
this upper bound includes many situations, it must take into account the case where both the 
velocity and the source vanish. Then c ~ ±1 as t — > oo, and by definition, ~ |^|^, which is in 



agreement with Eq. (4.14). These results enhance the existence theory obtained in Sec. 2.6, and 



give information about the long-time behaviour of solutions. 

As mentioned in Ch. [T| it is desirable in many applications to suppress concentration fluctuations, 
since this leads to a homogeneous mixture. In this report, we use advection as a suppression 
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mechanism, and we would therefore hke to know the maximum suppression achievable for a given 
flow. This inspires us to seek lower bounds on nip, in addition to the upper bounds found in this 
section. 



4.4 Lower bounds on the concentration fluctuations 



In this section we discuss the significance of the lower bound on the rrip measure of concentration 
fluctuations. Due to the powers of the concentration that appear in the Cahn-Hilliard equation, 
it is possible to obtain an explicit lower bound for m4, which we then use to discuss stirring 
as a mechanism to suppress concentration fluctuations. Using Holder's inequality, a flow that 
suppresses concentration fluctuations in the sense will also suppress them in the rrip sense, for 
p G [1,4], and thus our discussion of the suppression of concentration fluctuations by stirring is 
valid in this broader sense. 



As discussed in Sec. 5^, a suitable measure of concentration fluctuations for a symmetric mixture 
is 



rrir. 



0/ ' 



where c{x,t) is the concentration of the binary mixture and ||c||p is its norm. For p = 2, this 
gives the usual variance, used in the theory of miscible fluid mixing [29l ES, EQl |3T]. In that case, 
the choice p = 2 is a natural one suggested by the hnearity of the advection-diffusion equation. 
In the following analysis of the advective Cahn-Hilliard equation, it is possible to find an explicit 
lower bound for and we therefore use this quantity to study the suppression of concentration 
fluctuations due to the imposed velocity field. Given this formula, we can compare the suppression 
achieved by a given flow with the ideal level of suppression, and decide on the best strategy to 
homogenize the binary fluid. 



To estimate m^, we take Eq. (2.2), multiply it by an arbitrary, spatially-varying test function (x) 



and then integrate over space and time, which yields 



(4.16) 



where Q is the linear operator i> ■ V — -DA — ^DA"^. Using the constraint equation (4.16), the 
monotonicity of norms, and the Cauchy-Schwarz inequality, we obtain the following string of 
inequalities. 



< (l|c||2||Q0 + I)c^A0||2) < {\\cU)H\\Q<j) + Dc'A<f)r,)-^, 



[ 








Jn 





which gives the relation 



41)- > 



{WQcp + Dc^AcPWiy. 
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l)-^+D{\\c'Am)-^., 



We study the denominator 

(||g0 + Dc2A0||2)i < 

< 

where this bound follows from the triangle and Holder inequalities. Thus we have the result 



{Q<j>yd-xY +D\\A<p\U\\c\\t)-^.^ 



{J^{Q<Prd-x)'^ + D\\A<p\U\\c\\t) 



1 ) 

2 



or 



\c\\i) 



2\ i 

2 



{Q(f)yrx 











/ s(f)d'^x 


D\\A 


0l|oo(||c||^ 




> 










Jn 



Using the monotonicity of norms (4.4), we recast this inequality as one involving only a single 
power-mean, 



I^IMIIcll^)^ 



(Q^i^yrxy + D\\A(i)\u\\c\\t)-^ 



> 


/ s(f)d"-x 




Jn 



Therefore, we have the following inequality for m4 = (||c||4)4, 

1714 [qo{v,D,-^) + D\\A(j)\\^ml] > I^^P^ 

where 



schd^x 



Thus, we obtain a lower bound for the measure of concentration fluctuations 

m4 > mf^ D\\A<^\\^ (mf + go ^, 7) ^f" - 
The cubic equation satisfied by m™"^ has a unique positive root. 





d^'x 


Jn 





(4.17) 



To probe the asymptotic forms of (4.17), we rewrite the forcing terms v{x,t) and s (x) as an 
amplitude, multiplied by a dimensionless shape function. Thus, 



V = VoV , 

S = SqS , 



Vo = \n\-H\\v\\l)-2., 
So = ir^rziisiL. 



Then for a fixed value of 5*0 and D, and Vq ^ 1 (large stirring), we have qo ^ Vo{J {v ■ V0)^ d'^x) 2 , 



and the lower bound m™^° takes the form 
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On the other hand, for fixed 5*0 and Vq, and D ^ 1 (large diffusion), the lower bound takes the 
form 

mr~7^ — ^-'^ L , D^l. 

It is possible to obtain similar asymptotic expressions for m2, by minimizing a constrained func- 
tional of the concentration. Apart from a volume factor, the asymptotic form of m2 agrees exactly 
with the asymptotic form of just obtained. The functional to be minimized is 

$ [c] = y c^d^x^ ~Kj {'^^^ ^ Dc^^4> + s4>^ (Tx^ 

where (x) is a test function and A is the Lagrange multiplier for the constraint. This approach 
has been taken in [30] for the advection-diffusion equation. Setting 6^/6c = gives 



1 - \/l - 12A2L)A0Q0 



Evaluation of 6'^^/6c6c' shows that Eq. (4.18) produces a minimum of $ [c]. Given the expression 
Q = Vqv ■ V — -DA — ■yDA'^, the minimum Eq. (4.18) is XVov ■ V0 at large Vq. Substitution of this 
expression into the constraint 

' ) = 



cQ(f) + Dc^A(f) + s(, 

gives A = - (So/V^) [{J^ scpd^'x) / {J^ {v ■ V^f ci"x)] , and hence 

S'n I fo s(hd"'x I 

"^2 ~ TTTir— — — TT, V^O > 1- 

Vo {j^{v ■ Vc) rf"x)2 

For fixed Vq and 5*0 and large a similar calculation gives 

S'n I fo scbd^x I 
mT'" ~ ^ ^ r, D > 1 



D 



[j^ (A0 + 7A20)2 ^„^j 



1 ) 

2 



These expressions show that, apart from a volume factor, the lower bounds on the m2 and m4 
measures of concentration fiuctuations are identical in the physically interesting limits of high 
stirring strength, or high diffusion. In particular, the expression 



1 „:„ 5'o \j^s0d"-x\ 



mf", M-^mf'^-f ^ , ' for large K): 



indicates that if a flow can be found that saturates the lower bound m™™, the suppression of 
concentration fluctuations can be enhanced by a factor of Vq^ at large stirring amplitudes. Such 
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a flow would then be an efficient way of mixing the binary fluid. 



4.5 Scaling laws for 7724 



In this section we investigate the dependence of the 777,4 measure of concentration fluctuations on 
the parameters of the problem, namely the stirring velocity v, the source s, and the diffusion 
constant D. For simplicity, we shall restrict our interest to a certain class of flows, which enables 
us to compute long-time averages explicitly. 

The lower bound for the measure of concentration fluctuations is the unique positive root of 
the polynomial 



D\\A<P\\^{mf-y + qoiv,D,^)mf'^-\n\-^ / s^t/^x 

Jn 

where (j) (x) is a test function and 

qo {v, D,-f) = (^j [vV(p- DA(j) - D^A'^cf] ^ (Fx 



0, 



(4.19) 



We specialize to velocity fields whose time average has the following properties, 

(x, ■)) = 0, {vi{x,-)v.j{x,-)) 



n 



(4.20) 



The flow V {x,t) is defined on the 7i-torus [0,//]"". A statistically homogeneous and isotropic 
turbulent velocity field automatically satisfies the relations ( 4.20[ ), although it is not necessary for 
V to be of this type. The source we consider is monochromatic (that is, it contains contains a 
single spatial scale) and varies in a single direction. 



s = V2S0 sin {ksx) . 

Our choice of velocity field makes the evaluation of go (*^;-D,7) particularly easy; it is 



(4.21) 



go (^,^,7) 



1/2 

_iL 
n 



\V(f)frx + D^ / (A0 + 7A 



In studies of the advection-diffusion equation ^SOJ, it is possible to find an explicit test function (p 
that sharpens the lower bound on 77^2. The procedure for doing this depends on the linearity of 
the equation. Here, this is not possible, and for simplicity we set = s. This choice of certainly 
gives a lower bound for the 7774 measure of concentration fluctuations, with the added advantage 



of enabling explicit computations. Having specified the coefficients of the polynomial in Eq. (4.19) 



completely, we extract the positive root of this equation, and find the lower bound m^^^, as a 



function of Vq. The results of this procedure are shown in Fig. 4.1 



Let us examine briefly the scaling of the upper bound m^^^ with the problem parameters. The 
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Figure 4.1: (a) The lower bound for m™° as a function of Vq for a monochromatic source; (b) the 
dependence of m™™ on the velocity amplitude Vq. In (a), the scale of the source variation decreases 
in integer multiples from kg = 2tt/L in the uppermost curve, to kg = Sir/ L in the lowermost curve, 
while in (b) the source scale is set to 2tt/L. In both figures, we have set D = So = 1- 



upper bound satisfies the polynomial equation 



max ^4 



+ 



2\n\i 



D 



L 
2tt 



m 



max 



111 111 



\L J 



(4.22) 



which depends only on the diffusion D, the source amplitude 5*0, and the maximum stretching 

1 

Woo. For large 1^00, the flow-dependence of the upper bound is m™^"^ ~ W(^. This dependence is 
verified by obtaining the positive root of Eq. (4.22), which is as a function of Woo- The results are 



shown in Fig. 4.2 




Figure 4.2: The dependence of the upper bound m™^^ on the maximum stretching Woo- The source 
affects the upper bound only through its square mean Sq = \Q\~2\\s\\2. 
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Mixedness at small scales 



In the same manner, we find a lower bound for tlie quantity (|| VCII4) 4 . A similar gradient norm has 
arisen in the problem of quantifying mixedness at small scales in the theory of miscible fluids [30] • 



Using the constraint equation (4.16) again, we find 



x\ 



1 ; 

2 



(I|VC||^)5 > — , 

" (||Qv0 + 3c2DV0||i) 

where Qv = v — DV — D7AV. We study the denominator 

{\\Q^<P + 3Dc'vm)'^ < {\\Qvnl)^+SD{\\c'vm)K 

< (||Qv0||^)^+3Di|V0|U(||c^ll^)^ 



(4.23) 



{QycPfrxY +3D\\V(j)\U\\c\\t)l 

Using the Gagliardo-Nirenberg-Sobolev inequality on the periodic domain = [0, L]", we obtain 

n < 4. 

Thus, the denominator is bounded above by 



|c||4 < Cg||Vc||^ < Cg\^\ 2 ||Vc||2, 



{Q^^f<rxy + 3D||V0||oofco(||Vc||l)^ = clM"^-'. 



We therefore have the inequality 



|fi|MIIVc||^)i 



and hence 



|Vc||l)3>m, 3D\\V(j)\\ooknm^ + qo,s7iv,D,-f)m-\n\ ^ 



> 


I' 


(pcTx 




Jn 




/ 


bdTx 


= 0, 


Jn 







(4.24) 



where 



[v - DV(t) - L>7AV0]^ (Tx 



For a velocity field with the time-average relation {v-i (aj, ■) Vj (tc, ■)) = V^n ^6ij, this is 



go,v (v, 7) = <! — l^^l [ (V</) - 7AV<^)' 

^ Jn 



This concludes our investigation into lower bounds on the level of mixedness achievable in a 
stirred, phase-separating fluid. We note briefly that for any nonzero source, the lower bound m; 



mm 

■4 
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is nonzero, meaning that no matter how hard one stirs, there will always be some inhomogeneity 
in the fluid, and this is in fact true for any flow. However, the quantity m™™ tells us how much 
homogeneity we can achieve and is therefore a yardstick for stirring protocols. We use this yardstick 
to test model flows in the next section. 



4.6 Numerical simulations 



In this section we solve Eq. (2.2) numerically for two flows, and verify the bounds obtained in 



Sees. 4.3-4.5 We use the sinusoidal source term in Eq. ( |4.21 ) with periodic boundary conditions, 
and the source scale ks therefore takes the form (27r/L)j, where L is the box size and j is an 
integer. We speciahze to two dimensions and study two standard flows that are used in the 



analysis of mixing: the random-phase sine flow of Sec. 3.3, and the constant flow [31] . 



Random-phase sine flow 



The random-phase sine flow is the time-dependent two-dimensional flow 



t) = V2Vo sin {Ky + (pj] 



Vy {x, y, t) = V2V0 sin {k^x + ipj) , 



Vx 



JT<t< (j + I) T, 

(j + ^)r <t< (j + l)r. 



(4.25) 



where and ijjj are phases that are randomized once during each flow period r, and where the 
integer j labels the period. The flow is defined on the two-dimensional torus [0,L]^. The time 
average of this velocity field has the properties listed in Eq. (4.20). We solve Eq. ( 2.2[ ) with the flow 
in Eq. (4.25) as in Sec. 3.3 Recall that we used the lattice method for advection [231|27j, which 



was effective only when the spread of matter due to diffusion is much slower than the spread due to 
advection, that is, t/Tjj -C 1. Thus, as before, we therefore set D = 10~^, while we set r = L = 1. 
A numerical experiment with Vq = shows that 6*0 = 5 x lO""^ gives rise to a morphology that 
is qualitatively different from the sourceless morphology, and we therefore work with this source 
amplitude. Finally, following standard practice [161 ED], we choose 7 ~ Ax^, the gridsize. 

Using this nondimensionalization, and the identity Woo = V^Vok^ for the sine flow, we recall the 
large-stirring forms of m™'^^'™™. For Vq ^ the lower bound has the form 



So Jq sin (ksx) (pd'^x 



Vn 



(4.26) 



with the power-law relationship m™™ ~ ^, while for Vq ^ 1 the upper bound has the form 



2U„ 



1/2 
^^0 • 



(4.27) 
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These scaling results are identical to those for the advection-diffusion problem |28j . 



Before studying the case with flow, we integrate Eq. (2.2) without flow, to verify the effect of 



the source. For a sufficiently large source amplitude, the concentration phase separates and forms 







(a) 



(b) 



(c) 



(d) 




(e) 

Figure 4.3: The concentration field for Sq = 5 x 10"^ and (a) t = 100; (b) t = 1000; (c) t = 2000; 
(d) t = 8000. A steady state is reached in (d), evidenced by the time dependence of m4 in (e), 
where is constant for t ^ 1. 

domains rich in either binary fluid component. These domains are aligned with variations in the 
source. A steady state is reached and m4 attains a constant value, as seen in Fig. 4.3, On the 



other hand, for small source amplitudes, we have verified that that the domains do not align with 
the source, and their growth does not saturate. We do not consider this case here, since we are 
interested in sources that qualitatively alter the phase separation. These different regimes are 
discussed in 



We consider the case with flow by varying Vq, and find results that are similar to those found in 



Sec. |3.5[ for the same stirring mechanism without sources. For all values of Vq, the concentration 
reaches a steady state, in which ||c||4, effectively the pre-averaged form of m4, fluctuates around 
a constant value. For small values of Vq, the domain growth is arrested due to a balance between 
the advection and phase-separation terms in the equation, while for moderate values of Vq, the 
domains are broken up and a mixed state is obtained. At large values of Vq, the m4 measure of 
concentration fluctuations saturates: further increases in Vq do not produce further decreases in 
m^. At these large values of Vq, the source structure is visible in snapshots of the concentration, 
as evidenced by Fig. |4.4[ 



We investigate the dependence of concentration fluctuations on the stirring strength Vq, and show 
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the results in Fig. 4.5 The theoretical upper and lower bounds on depend on Vq and are 



obtained as as roots of Eqs. (4.19) and (4.22). In the limit of large Vq, these bounds have power- 



law behaviour, with m™*^^ 



^0 



and m™™ 



as demonstrated by Eqs. (4.26) and (4.27). The 



numerical values of m4 are indeed bounded by these limiting values, although the Vo-dependence 
is not a power law. Instead, the function (Vq) is a nonincreasing function, with a sharp drop 





i:: 




\^ 10 



(b) 



(c) 



Figure 4.4: A snapshot of the steady-state concentration field for (a) Vq = 0.001; (b) Vq = 0.1; 
(c) Vq = 10. The source parameters are Sq = 5 x 10~^ and kg = 27t. In (a) domain growth is 
arrested, in (b) the domains are destroyed and the binary fluid mixes, while in (c) m,^ measure of 
concentration fluctuations is minimized, and the source structure is visible. 



occurring in a small range of Vq- values. Thus, the fluid becomes more homogeneous with increasing 
Vq. We discuss the effect of stirring on the inhomogeneity of the fluid by introducing the notion 
of mixing efficacy. 

We measure the ability of a given stirring protocol to suppress concentration fluctuations by the 
mixing efficacy. Similar ideas are often applied to the advection-diffusion equation 
We define the dimensionless mixing efficacy 



mf-^ {Vo = 0) 



iVo) 



For a given flow, the number rjp quantifies the flow's ability to suppress concentration fluctuations. 
In a well-mixed flow, the local deviation of c{x,t) away from the mean will be small; a small 
-value is a signature of a well-mixed flow. We are therefore justified in calling rip the mixing 



efficacy. We obtain some control over the mixing efficacy r]4 from the inequalities of Sees. |4.2 
and |4.4[ Based on these inequalities, the mixing efficacy is bounded above and below. 



^min _ "^4 
'14 



{Vo = 0) 



(Vo. 



<V4< 



mf ° {Vo = 0) 



(Vo) 



max 

'k ■ 



We have plotted the upper and lower bounds on the mixing efficacy for the case of monochromatic 



sources in Fig. 4.5 (b). The maximum efficacy always exceeds unity in this case, which implies 
the possibility of finding stirring protocols that homogenize the fluid. On the other hand, the 
minimum efficacy is less than unity, which indicates the possibility of finding stirring protocols 
that actually enhance concentration fluctuations, and this enhancement depends weakly on the 
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Figure 4.5: (a) The m,^ measure of concentration fluctuations or mixing for the sine flow, as a 
function of the stirring parameter Vq. The values of D and 5*0 are given in the text. The upper 
and lower bounds are shown for comparison; (b) The mixing efficacy for the sine flow, with the 
upper and lower bounds shown for comparison. 



maximum stretching Woo- This latter case is not surprising, given that a uniform shear flow causes 
the domains of the Cahn-Hilliard fluid to align, rather than to break up. The sine-flow efficacy is 
a nondecreasing function of the stirring parameter Vq. At small values of Vq, small increases in the 
vigour of stirring lead to small increases in the mixing efficacy. There is a window of intermediate 
Vq- values for which the mixing efficacy increases sharply with increasing Vq. At higher values of 
Vq, the efficiency saturates, so that further increases in the vigour of stirring have no effect on 
concentration fluctuations. In the saturated case, the concentration field mirrors the structure of 
the source function, as in Fig. 4.4 (c). 



Constant flow 



We study the flow {vx,Vy) = (Vo,0), where Vq is a constant. We choose a nondimensionalization 
that is set by the diffusion time Td = L^/D, and obtain the following parametric version of 



Eq. (2.2), 



where Vq 



LVo/D, r 



dc ,dc 
dV ^dx' 

and 5'( 



A' {c^ -c- 7'A'c) + V2S'q sin {k[x') 



(4.28) 



S'qL'^/D. We immediately drop the primes from 
Eq. (4.28). We fix 7 and 5*0 and vary the flow strength Vq. As in the case of no flow, we choose 



5*0 such that the morpohology is qualitatively different from the sourceless one; in the units used 
here, we set Sq = 50. This flow does not satisfy the time-correlation relations (4.20), although the 
maximum stretching has the simple form Woo = 0. The upper bound obtained in Eq. (4.14) is 
therefore independent of the flow strength. 



We solve Eq. (4.28) numerically for various values of Vq and present the results in Fig. 4.6 For 



small Vq, the morphology of the concentration fleld is similar to the flowless case seen in Fig. 4.3 (d). 



except now the domains are uniformly advected in a direction perpendicular to the source variation. 
The small- Vo case is shown in Fig. 4.6 (a). As Vq increases, the domain boundaries are distorted 
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(a) 



(b) 



(c) 



Figure 4.6: A snapshot of the steady-state concentration field for (a) Vq = 10; (b) Vq 
Vq = 1000. The source parameters are So = 50 and ks = 27t. 



100; 



due to the advection, while for large Vq, the advection is sufficiently strong to precipitate the 
formation of narrower lamellae, although the domain structure persists. 

The m4 measure of mixedness extracted from the numerical simulations is almost constant across 
the range of stirring parameters Vq G [0, 1000], and resides in the range defined by the bounds in 
Sec. |4.4[ For Vq = 1000, m4 is slightly smaller than its value at Vb = 0, due to the presence of more 



interfaces. This difference is small however, and increasing Vq does little to mix the fluid. This is 
not surprising, given that local shears are necessary to break up domain structures [20], and that 
such shears are absent from constant flows. What this example shows however, is the difference 
between a miscible mixture, and a phase-separating mixture. For a diffusive mixture with the 
sinusoidal source we have studied, the constant flow discussed here is optimal for mixing |3T]; for 
a phase-separating mixture, the constant flow badly fails to homogenize the mixture. 



4.7 Summary 

We have introduced the advective Cahn-Hilliard equation with a mean-zero driving term as a 
way of describing a stirred, phase-separating fluid, in the presence of sources and sinks. By 
specializing to symmetric mixtures, we have studied a more tractable problem, although one with 
many applications. 

Our goal was to investigate stirring protocols numerically and analytically, and to determine the 
best way to break up the domains in the Cahn-Hilliard fluid and achieve homogenization. To 
this end, we introduced the mp measure of concentration fluctuations. Since in a well-mixed fluid, 
the concentration exhibits small spatial fluctuations about the mean, with better mixing leading 
to smaller fluctuations, we used mp as a measure of mixedness or homogeneity. We proved the 
existence of mp for long times, for p G [1, 4], and obtained a priori upper and lower bounds on m^, 
as an explicit function of the imposed flow v {x, t), and the source s {x). 

We compared the level homogeneity achieved by the random-phase sine flow and the constant 
flow with the lower bound, and found that the sine flow is effective at homogenizing the binary 
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fluid, while the constant flow fails in this task. This is not surprising, since it is known that 
differential shears are needed to break up binary fluid domains, although it is radically different 
from the advection-diffusion case, where the constant flow was the optimal mixer. The question 
of whether or not a flow is a good mixer in this context was discussed using the mixing efficacy, 
defined in Sec. 4.6[ Given such a definition, it is possible to compare stirring protocols and find the 
optimal protocol for mixing a binary fluid. Our upper bound on the efficacy provides a meaningful 
notion of this optimality. This result may be useful in applications where the homogenization of 
a binary fluid is desirable, since we have set a lower limit on precisely how much homogeneity can 
be achieved. 



Chapter 5 



Nonlinear dynamics of phase separation 
in ultra-thin films: Analysis 



5.1 Overview 



In this chapter we recall the Navier-Stokes Cahn-Hilliard (NSCH) equations derived in Sec. 2.3 



shifting focus from the passive tracer, to the active one. We study these equations in the incom- 
pressible limit. We shall specialize to thin films, in which the binary liquid forms a thin layer on a 
substrate. Then, if the vertical gradients are small compared to horizontal ones, a long-wavelength 
version of the NSCH equations can be obtained, representing a much simplified problem. We 
present the derivation of these long-wavelength or thin-film Stokes Cahn-Hilliard equations. We 
analyze these equations and obtain results concerning the existence, regularity, and uniqueness of 
solutions. 



5.2 The thin-film Stokes Cahn— HiUiard equations 



In this section we recall the incompressible Navier-Stokes Cahn-Hilliard equations from Sec. 2.3 
and discuss the passage from these equations to the long-wavelength approximation in two dimen- 
sions. We enumerate the scaling rules necessary to obtain the simplified equations. Finally, we 
arrive at a set of equations that describe phase separation in a thin film acted on by arbitrary 
body forces. 



The incompressible NSCH equations derived in Sec. |2.3| describe the coupled effects of phase 
separation and flow in a density-matched binary fluid. If v is the fluid velocity and c is the 
concentration of the mixture, where c = ±1 indicates total segregation, then these fields evolve in 
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the following way, 

dv 1 

— + vVv = V-T V0, (5.1a) 

dt Po 

dc 

— + vVc = DA(c^-c--fAc), (5.1b) 
ot ^ ' 

V-v = 0, (5.1c) 

where 

T- ^(5- •+ { +^^"'^ /3 (5 2) 

is the stress tensor, p is the fluid pressure, is the body potential and po is the constant mixture 
density. The constant v is the kinematic viscosity, v = r]/po, where t] is the viscosity. Additionally, 
/? = eo/po is a constant with units of [Energy] [Mass] ^/^ is a constant that gives the typical 
width of interdomain transitions, and D a diffusion coefficient with dimensions [Length]^ [Time] 

If the system has a free surface in the vertical or 2;-direction and has infinite or periodic boundary 
conditions (BCs) in the lateral or x-direction, then the vertical BCs we impose are 

u = w = Cz = Czzz on z = 0, (5.3a) 

while on the free surface z = h{x,t) they are 

dT 

fiifijTij = Fk, hitjTij = (5.3b) 

dh dh o N 

^=-+.-, (5.3c) 

hidiC = 0, hidiAc = 0, (5.3d) 

where n is the unit normal to the surface, t, is the unit vector tangent to the surface, s is the 
surface coordinate, F is the surface tension and k is the mean curvature, 

K = V • n 

[Ox, Oz 



This choice of BCs guarantees the conservation of the total mass and volume, respectively 



M= / c{x,t)d'x, V= / d'x. (5.4) 

JBomit) JDomit) 

Here Dom [t) represents the time-varying domain of integration, whose time dependence is due to 



the variable nature of the free surface height. Note that in view of the concentration BC (5.3d) 



the stress BC (5.3b) does not contain c{x,t) or its derivatives. 
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Figure 5.1: Definition sketch for tlie long-wavelengtfi approximation in two dimensions. Tlie 
approximation relies on tlie smallness of the parameter 6 = ho/X. 



These equations admit a great simplification if the fluid forms a thin layer of mean thickness ho, for 
then the scale of lateral variations A is large compared with the scale of vertical variations h^. This 
is shown in Fig. 5.1 Specifically, the parameter 6 = ho/X is small, and after nondimensionalization 
of Eq. (5.1 ), we expand its solution in terms of this parameter, keeping only the lowest-order terms. 

For simplicity, we work in two dimensions. 



For a review of this method and its applications, see |50] 
but the generalization to three dimensions is easily effected. 

In terms of the small parameter S, the equations nondimensionalize as follows. The diffusion 
timescale is to = X^/D = h"^/ (S^D) and we choose this to be the unit of time. Then the unit 
of horizontal velocity is uq = A/to = SD/ho so that u = {6D/ho) U, where variables in upper 
case denote dimensionless quantities. Similarly the vertical velocity is w = {6'^D/ho) W. For the 
equations of motion to be half-Poiseuille at 0(1) (in the absence of the backreaction) we choose 
p = {riD/h'^)P and (f) = {riD/h'^)(^. Using these scaling rules, the dimensionless momentum 
equations are 



5Re 



d 



dX 
^uD dX 



(P + $) + 6'- 



+ 



dc 
dX 



+ 



dX^ 
dc 
dZ 



2 1 



P'y dc 
l^dX 



+ 



dX^ dZ 



(5.5) 



fdW ^dW BW 
I -t::^ + U—— + W 



dX 



d 



dZ J dZ 

I Pi d 



(P + $) + + 5^ 



2z/P) dZ 



\dX 



dX^ 
f dc 



dZ^ 

(3^ dc 
^dZ 



d'^c d'^c ' 



(9X2 



(5.6) 
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where 

fle=^ = «^ = 0(l). (5.7) 

The choice of ordering for the Reynolds number Re allows us to recover half-Poiseuille flow at 
O (1). We delay choosing the ordering of the dimensionless group jS'-f/Du until we have examined 
the concentration equation, which in nondimensional form is 

dc dc dc \ 

where = 'y/h^. By switching off the backreaction in the momentum equations (corresponding 
to Pj/Du — > 0), we find the trivial solution to the momentum equations, U = W = dx {P + ^) = 
dz {P + ^) = 0, H = 1. The concentration boundary conditions are then cz = czzz = on 
Z = 0,1 which forces = so that the Cahn-Hilliard equation is simply 

^ = ^(c'-c)-s'r'— 

To make the lubrication approximation consistent, we take 

5C„ = a = W^o = 0(l). (5.9) 



We now carry out a long-wavelength approximation to Eq. ( 5.8[ ), writing U = Uq + O {5), W 



Wq + O (5), c = Cq + 5ci + S^Ci + .... We examine the boundary conditions on c (a;, t) first. They 
are n ■ Vc = n ■ VAc = on Z = 0, iJ; on Z = these conditions are simply dzc = dzzzc = 0, 
while on Z = if the surface derivatives are determined by the relations 

n • V oc -6^Hxdx + dz, 

h ■ VA oc -5'^Hxdxxx - S^Hxdxdzz + 5'^dxxdz + dzzz- 

Thus, the BCs on cq are simply dz(^ = dzzzCo = on Z = 0, H, which forces cq = co(X, T). 
Similarly, we find Ci = Ci (X, T), and 

dc2 ryHxdco &^C2 HxdcQ 7 ^ m cri 

^ = ^irdX^ dZ^ = irdX' ^or^nyZe[0,H]. 



In the same manner, we derive the results dzzzzC2 = dzzzzCs = 0. Using these facts, equation (5.8 ) 
becomes 



dco dco 



gX2 (^0 Co) Cl^^^ + {3cl 1) ^ ^^'dx^ H dX ^^dZ^- 
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We now integrate this equation from Z = to Z = H and use the boundary conditions 



dz^ 



on Z = 0, 



d^co d f Hxdco 



H- 



dz^ " ' ^dx \H dxj dx^\H ax 

After rearrangement, the concentration equation becomes 



on Z = H. 



^df^^^'dX 



H 



dX^ 



+ 



OH d 
dXdX 



Co - Q 



^2 ^^Cq ~ 2 Hx dct 



a 



dX^ ^ H dX 



where 



if 



H 



Uo{X, Z,T)dZ, 



is the vertically-averaged velocity. Introducing 



= Co - Co - 



the thin-film Cahn-Hilliard equation becomes 



C' d 



H dX \ dX 



H 



dco 



OT ^dX ~ Hdx[ dX 



(5.10) 



We are now able to perform the long- wavelength approximation to Eqs. (5.5) and (5.6). At lowest 



order, Eq. (5.6) is 



d 



dco 



(P + + r- 
dZ^ ' dz 



1 d f dcp 

^""dz [dx 



0, 



which is simply dz (P + $) = 0, and hence 

P + $ = Psurf + $surf = P {X , H (X, T) , T) + $ (X, H (X, T) , T) . 

Here 



Dp 



Du 



0(1] 



is the dimensionless measure of the backreaction strength. Thus, Eq. (5.5) becomes 



dZ^ ~ dX ^ ^ ^"'dXKdX ^''dxdz^-^ 



(5.11) 



5.2. The thin-Glm Stokes Cahn-Hilliard equations 



94 



and using dzzC2 = (Hx/H) (dco/dX) this is 



d^Uo d 



r d 

^(Psurf + $surf) + ^^ 



H 



dco 
dX 



(5.12) 



At lowest order, the BC (5.3b) becomes 

dUo _ dt 
~dZ ~ dX 



on Z = H, 



(5.13) 



where T is the dimensionless, spatially- varying surface tension. Combining Eqs. (5.12) and (5.13) 
yields the relation 



^ = ^ + (^-//)^(Ps.. + $s..) + ^^ 



H 



dco 
dX 



Making use of the BC Uq = on Z = and integrating again, we obtain the result 



(5.14) 



The vertically-averaged velocity is therefore 



H 



dco 
dX 



Using the standard Laplace- Young free-surface boundary condition, this becomes 



1 d^H 



\dX\ CdX^ 



surf 



r d 
HdX 



H 



dco 
dX 



where 



C 



upD 



0(1) 



(5.15) 



/ioo-o52 

Finally, by integrating the continuity equation in the Z-direction, we obtain, in a standard way, 
an equation for free-surface variations. 



Let us assemble our results, restoring the lower-case fonts and omitting ornamentation over the 
nondimensional quantities, 

dh dJ 
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dx 



dx 



where 



J 



dx ^ \dx 



1 d^h 
C dx' 



+ 



+ 



r d 
h dx 



h(—\ 
\dx J 



= c^-c-^ — 
^ h dx 



^ dx) 



(5.16b) 



(5.16c) 



(5.16d) 



and where we have the following nondimensional constants 



c 



upD 



Dv ho hoaoS"^ 

These are the thin-film Stokes Cahn-Hilliard (SCH) equations. The integral quantities M and V 



defined in Eq. (5.4) are manifestly conserved, while the free surface and concentration are coupled. 



We note that the relation = ^y/l /h^ = O (1) is the condition that the mean thickness of the 
film be much smaller than the transition layer thickness. In experiments involving the smallest 
film thicknesses attainable (10~^ m) |lUlj . this condition is automatically satisfied. The condition 
is also realized in ordinary thin films when external effects such as the air-fluid and fluid-substrate 
interactions do not prefer one binary fluid component or another. In this case, the vertical extent 
of the domains becomes comparable to the film thickness at late times, the thin film behaves in a 
quasi two-dimensional way, and the model equations are applicable. 

The choice of potential determines the behaviour of solutions. If interactions between the fluid 
and the substrate and air interfaces are important, the potential should take account of the Van 
der Waals forces present. A simple model potential is thus 

= Ah-P, 

where A is the dimensionless Hamakar constant and typically p = 3 [50]. Here A can be positive 
or negative, with positivity indicating a net attraction between the fluid and the substrate and 
negativity indicating a net repulsion. This choice of potential can also have a regularizing effect. 



preventing a singularity or rupture from occurring in Eq. (1.1 ). 



For (f) = —\A\/h^ (repulsive Van der Waals interaction), the system of equations (1.1) has a 
Lyapunov functional F = Fi + F2, where 



dx 



J^fdhV \A[ 
2C \dx) 2 1^2 



r 

CI 



dx h 



dc 
dx 



and where Q = [0, L] defines the lateral extent of the system. In this chapter we take Q to 
be a periodic domain, which fixes the boundary conditions, although other choices of boundary 
conditions are possible. By differentiating these expressions with respect to time, we obtain the 
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relation 



F1 + F2 



3 ^^(ix/i '^g^ 1^(7(9x2 ^ J hdx 



h{-\ 
dx J 



dxh . 



(5.17) 



which is nonpositive for nonnegative h. This fact holds the key to the analytic results obtained in 
the next section. 



5.3 Existence and regularity of solutions 



In this section we prove that solutions to the model equations do indeed exist. We focus on the 
one-dimensional equation set 



dh 

di 



d_ 
dx 



f{h) 



dx'^ 



+ 



dx 



1 dh 
g (h) dx 



, d ffjh) d 
dx\ g (h) dx 



d_ 

Ft 



(cgm 



where 



d_ 
dx 



cf{h) 



d^h 
dx^ 



+ 



d_ 
dx 



c dh 



^<^[J{h) d 



g{h)dx\ dx\ g (h) dx 

-—Igih) — 

dx \ dx 



9{h) 



1 d 



dc 
dx 



c — c — 



g{h)dx V^^^^<9x) 



f{h) = h\ g{h) = h. 



(5.18) 



The application we have in mind involves a periodic domain Vt = [0,L]. We shall prove our result 
for smooth initial conditions, 



h{x,0) 



hoix) > 0, c(x,0) = Co(x), (/io(x),Co(a;)) G H''' (n) . (5.19) 



We shall prove that the solutions to this equation pair exist in the strong sense; however, we shall 
need the following definition of weak solutions: 



A pair {h, c) is a weak solution of Eq. (5.18) if the following integral relations hold, 

-To 



dt / dx ipth = 
Jn 

dt / dx^J -/(/i) — + ^TT— + ^ ' 
Jn I 



dx^ g (h) dx g (h) dx 



gih)—(-\ 

dx \dx J 
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/ dt dx i/jfcg (h) = 
'0 Jn 

dt I dxijj -cfih)^ + + ^ 
'0 Jq I 

"To 



dx^ g (h) dx g (h) dx 



9{h) 



d [ dc 



dx \ dx 



+ I dt J dxij^^g{h) ^ 



3 1 d f dc\ 



(5.20) 



where Tq > is any time, and ip{x,t) and ilj{x,t) are arbitrary differentiate test 
functions that are periodic on Q and vanish at t = and t = Tq. 

We address the following statement, 



Given the initial data in Eq. (5.19), the equations (5.18) possess a unique strong solution 



endowed with the following regularity properties: 

{h, c) e (0, To, (Q)) n (0, To; H^'^ (Q)) n Ci'i {Q x [0, To]) . 



Recall that the function / resides in the class H'^'^ [Q) if 



.July <oo. 



The proof is given in the following steps. 



5.3.1 Regularization of the problem 



We introduce regularized functions fs (s) and gs (s) such that 



hm /,(.) = /(.) 

e— >0 



lim 0^(5) =g{s). 

£^0,S>0 



For now we do not specify fe (s), although we mention that a suitable choice of fs (s) will cure the 
degeneracy of the fourth-order term in the height equation. On the other hand, we require that 
the function g^ (s) have the following properties: 

• gs (s) = s + e, for s > 0, 

• gs (s) > 0, for s < 0, 

• lim^^.oo gs (s) = \e, 



gs (s) has as many derivatives as necessary. 
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Figure 5.2: A sketch of the regularized function (s) obtained in Eq. (5.21) 



One way of obtaining such a function is to define 



e 
2' 



s < — e, 



e + s, 



s > 0, 



(5.21) 



which is in the variable s, and has Lipschitz third derivative; inspection of Eq. (|5.18l) shows that 



this degree of smoothness is sufficient to regularize the equations. A sketch of this regularization 



is shown in Fig. 5.2 



The regularized PDEs we study are thus 



ht 



Je = fe {h) K 



1 , fe (h) , 2^ 
-h^ — [ge{h) c^)^ 



where 



- - (cJe)^ - {9e (h) fle,x)^ , 

3 1 

He = C - C 

9e 

The c-equation can also be written as 

ct ■■ 



9e{h) 



9e (h) 9e (h) 



{9e (h) C^.)x • 



(5.22) 



JeCx H {9e {h) He 



9e (h) 9e (h) 



9s (h) 



[Je,x + 9e (h) ht] 



5.3.2 The Galerkin approximation 

We choose a complete orthonormal basis on the interval Q, with periodic boundary conditions. Let 
us denote the basis by {0j (x)}jgNo- We consider the finite-dimensional vector space span{0o, ■■■(Pn}- 
For convenience, let us take the (pi (x)'s to be the eigenfunctions of the Laplacian on [0, L] with pe- 
riodic boundary conditions, with corresponding eigenvalues — A^. Moreover, let 0o be the constant 
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eigenfunction. We construct approximate solutions to the PDEs (5.22) as finite sums, 

N N 
hN {X, t) = ^ VN,i (t) (pi (x) , Cn {x, ^) = ^ lN,i (t) (x) . 

i=0 i=0 

If tlie (smootli) initial data are given as 

oo oo 

h (x, 0) = /lo {x) = ^ r7°0i (x) > 0, c (x, 0) = cq (x) = 7°0i (x) , 

i=0 1=0 

then the initial data for the Galerkin approximation are 

N N 

hN {x, 0) = h% (x) = T]i(j)i (x) , Cn (x, 0) = c% (x) = 7°0i (x) , 

i=0 1=0 

and the initial data of the Galerkin approximation converge strongly in the {Q) norm to the 
initial data of the unapproximated problem. Thus, there is a A^o ^ such that h% (x) > 0, 
everywhere in Q, for all > A^o- Henceforth we work with Galerkin approximations with > A"o. 

Substitution oi Hn = Yl!i=QVN,i4>i into a weak form of the /i-equation yields 

|((/i^,</),)) = ((4^,0,g), (5.23) 
where ^ 

Je {hN, Cn) = fe (hN) hN,xxx 77-^hN,x ^77^ ids (^n) 0% ) , 

is the flux for the regularized /i-equation, and in this chapter, ((0, ip)) denotes the pairing (jfijjdx. 
This is recast as 

where the function $jv {vn,1n) depends on 

Vn = {vn,o, ■■■,riN,N) , In = {1n,o, ■■■,1n,n) , 

and is locally Lipschitz in its variables. This Lipschitz property arises from the fact that the 
regularized flux, evaluated at the Galerkin approximation, is a composition of Lipschitz continuous 
functions, and therefore, is itself Lipschitz continuous. 

Similarly, substitution of cat = J^iLo lN,i4'i into the weak form of the c-equation yields 

^ {{g Qin) Cn, (pj)) = {{Ke,N, , (5.24) 
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where 



Ke {hN, Cn) = CNfe (/^Af) hN,2 



cn 



9e {h 



-h 



N 



9e [hN) 



is the flux for the regularized c-equation. Rearranging gives 

{{9 {h^) cn,u <Pj)) = {{Ke.N, 4>j,x)) - {{9' {hn) CNhN,t, 4>j)) , 
and the left-hand side can be recast in matrix form as Yla=o^oij'iN,a, where 



M 



VN,i(Pi <P 



'31 



which is manifestly symmetric. It is positive deflnite because given a nonzero vector (.^o, ■■■^iN)-, 
we have the relation 



a,j \ i J aj 



> for {^o,...,^n)^0, 

which follows from the positivity of the regularized function (s). We therefore have the following 
equation for •jM^i (t), 



dt 



N 



J2 {{K,^N, <Pa,x)) - {{9's (hN) CnhN,U 0a)) 



(5.25) 



a=0 



Inspecting the expression for Maj, 9 {hj^), and K^^^, we see that rj^- and ■jn' appear in a (locally) 
Lipschitz-continuous way in the expression J2a=o^al {{^e,N,<Pa,x)), while owing to the imposed 
smoothness of (s), the variables rj^, 77V and fjiy = {i)n,o, ■■■,Vn,n) appear in a Lipschitz-continuous 
way in the quantity 

AT 

The vector fj^ can be replaced by the function {vn,1n) and thus we obtain a relation 



in place of Eq. (5.25), where ^N,j (vn^In) is Lipschitz. We therefore have a system of Lipschitz- 
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continuous equations 

= ^N,j iVN, In) , = ^7v,i iVN, 7iv) , 

and thus local existence theory guarantees a solution for the rjj^i/s and 7Ar,i's for all times t in a 
finite interval < t < a. This solution is, moreover, unique and continuous. 



5.3.3 A priori bounds on the Galerkin approximation 

We identify the free energy 

F = [ dx [Ihl + (h)] + [ dxg, (h) h {c' - 1)' + Icl] , 
where G" (s) = 1/ [fs [s) (s)]. Since the Galerkin approximation satisfies the weak form of the 



PDEs given in Eqs. (5.23) and (5.24), it is possible to obtain the following free-energy decay law, 
dF 



dt 



dx fe {h 



N 



-h 



N,xxx 



+ 



h 



N,x 



+ 



fe (hN) 9e (hN) 9e {hN 



{9e (h^) C% ,) 



1 2 



- / dxge^h^) I dxJe,N 

'n Jn^ 



-h 



N,xxx 



+ 



h 



N,x 



fe ihN)9e (/^TV). 

+ J dx^^g'^ {h]y) J^^N (c^ - C - Cn,xx) - ^■e,NCN [q's {hN) + Je,N,x] " Je,N ^-e^NCN , 

< t < a, (5.26) 

where f2_ {t) = {x G ^7|/^A^ (x, t) < 0}. Now given the time continuity oih^ (x, t) in (0, a), and the 
initial condition /i^ (x) > (since > Nq), there is a time ai, such that < cti < a, and such 
that /iTv {x, t) > for all x G and all t G (0, Therefore, (t) = for t G (0, ai), and hence 

F [cn (x, t) , Hm (x, t)] < F [cn (x, 0) , (x, 0)] < sup F [cm (x, 0) , hf^i (x, 0)] < oo, 

£,Afe[0,oo) 

for < t < (Ti. Consequently, we obtain the bound ||/iAr,a;||2 < ^i, where < t < cti, and where ki 
depends only on the initial conditions. We have Poincare's inequahty for /iat^x, 



\hN\\l< 



dxh^ [x 



+ 



L \ 



2tx I 



1 II^Af,zl|2- 



Now dx (x, t) = LrjNfl if)- Inspection of Eq. (5.23) shows that rjNfl if) = 1]^^ (0) = rf^. Thus, 



hNWi < L'W+ ( — J h = k2< oo. 
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Using the following inequality from Sec. 2.4 



1 1 „ „ ^„ „ 



we obtain the bound ^ 

||^Af||cxD < ^^ll^jvlb + V^||/lAr,x-||2, 

VL 

and hence ||/iAr||oo < ^3- Additionally, the following properties hold: 
• The function hj\f is Holder continuous in space, with exponent |, 

where these results hold in < t < ui, and where the constants ki, ^2, ^3, and k^ are independent 
of e, N, a, and ai, and in fact depend only on the functions (x) and Cq (x). 

In order to continue the estimates to the whole interval (0, a), we need to prove that hj\f (■, > 
almost everywhere (a.e.). If this is true, there is a new interval [(Ji,(T2), cti < < cr, on which 
Hn {-jt) > a.e., and we can then provide a priori bounds on Hn {-jt) and cat on the interval 
[ci, 0-2). It is then possible to show that (■•, 0^2) > a.e. and thus, by iteration, we extend the 
proof to the whole interval (0, a), and find that (., t) > a.e. on (0, a). 

We have the bound 

dxGeihN{-,t)) <ki, (5.27) 



a 

where ^4 depends only on the initial conditions, and where < t < ai. We now specify G^ (s), in 
more detail. This function satisfies the condition 

fe [S) Qe [s) 

We take (s) to be as defined previously, while fe {s) can be regularized as Qe {sf , which is 
Lipschitz continuous. By defining 

/"°° dr 

Ge{s) = - GM = - drGeir), 

Js fe [r) ge [r) Js 

we obtain a function G^ (s) that is positive for all s G {—00, 00), and 

Using the boundedness of G^ (s), and the time continuity of /i„ (■,)!:), we employ the Dominated 
Convergence Theorem, 

lim / (ix Ge (/lAT (■, t)) = / dx lim (/itv (■, ^)) = / dx G^ {h^ cri)) ^ ki. 
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Similarly, since the constant ki in the inequality ||/iAr,x||2 < fci, < t < (Xi depends only on the 
initial data, we extend this last inequality to t = (Xi, and thus hjsf {x,ai) is Holder continuous in 
space. 

In the worst-case scenario, the time ai is the first time at which [x, t) touches down to zero, and 

and we therefore assume for contradiction that {xq, cti) = 0, and that {x, cti) > elsewhere. 

1 

Then by the Holder continuity, for any x E Q we have the bound < (a^, cti) < ki\x — xq\^ , 
and thus 



dx Ge {JlN (■, (Ti)) 



dx 



n [hN {x, ai) + ey 



> 



6 



dx 



|x — Xol 2 + e 



^ - 6 



dx 



\x — xq\ + e [ 2y'L + e 



Hence 
6 

h 



dxGe {hN (^CTi)) 



> -2 loe 



£(2\^ + £^ +log| L-XO+ (2y/L + e^€ xq + (2^/1 + e | 



Thus, the integral J^Gs {hN {x,ai)) dx can be made arbitrarily large, which contradicts the e- 
independent bound for this quantity, obtained in Eq. (5.27). We therefore have the strong condition 



that set on which /^at (■, < is empty. Iterating the argument, we have the following important 
result: 



The set on which hN {-yt) < is empty, for < t < cr. 

Using the same argument, we have an estimate on the minimum value of hN {x,t), 

hr, 



oram^ Uliu hN{x,t) 

xen,te(o,cr] 



namely, 



^min + £ > -kiV L + \ kfL + 



kfL 



for all small positive e. Thus, 



hmin > B := -kiV L + \ kfL + 



kfL 



{5.21 



a result that depends only on the initial data co{x) and hQ{x). Now, using Eq. (5.28) and the 
boundedness result 



dx {hN) \ {c% - l) + 



2 , 1.2 



< k. 
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where depends only on the initial data, we obtain a priori bounds on \\cn. 



2 



It is also possible to derive a bound on ||cjv||2- We have the relation 



/ dx\{cl-lY< 
Jn 



which gives the inequality 1 1 Cat I II < 2||cAr||2+(4/i;5/-B) . Using the Holder relation ||cAr||2 < 4 ||cAr||4, 
we obtain a quadratic inequality in HcatHI, 



\cN\\i<2m\\c^\\i + ^^ 



with solution 

\\cn\\2 < M H ^> 

as required. From the boundedness of ||cAr^x||2 and ||civ||2 follows the relation ||cAr||oo < fcg < cc, a 
result that depends only on the initial conditions. Let us recapitulate these results: 

• II^Af.i'lh, is uniformly bounded; 

• 1 1 ^Af 1 1 00 is uniformly bounded; 

• The function \s Holder continuous in space, with exponent |; 

• The function h]^ is nonzero everywhere and never decreases below a certain value i? > 0, 
independent of A^, e, and a. 

• ||cAr,a;||2 is Uniformly bounded; 

• ||cAf||oo is uniformly bounded; 

• The function cat is Holder continuous in space, with exponent |, 
where these results are independent of A^, e and cr, and hold for < t < cr. 

5.3.4 Equicontinuity of the Galerkin approximation; convergence of 
Galerkin approximation 



Using Eq. (5.26), we obtain the bound 



dt' I dx\fe {Hn) 





h]^ X 1 / /7^2^ 

-nN,xxx H /, w /, N H 77—: [9e [hN) CATxj, 

9s{hN) Je{hN) ge{hN) 



-\ 2 



<F{0), 0<t<a, 
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a bound that is independent of N, a, and e. Since the quantity (/liv) ||oo = (||/iAr||oo + 
is bounded above by a constant Ai independent of N, e, and a, we have the following string of 
inequalities, 



dt' / dx J^ j^ 
Jn 



dt' / dx\fe {Hn) 
Jn 



-hN,xxx H /, \'f fi, \ TTTT [9e [hN) Cn,x)^ 

9e{hN)je{hn) ge{K) ' 



< / dt'\\fe{hN)\\oo / dx{h{hN) 

Jo Jn I 

<Ai f dt' I dxifeih^) 

Jo Jn I 



^N,xxx ~l~ 



hN,: 



+ 



9e {h^) fe {hn) ge 
tlN.x 1 



(7^ ids (hN) C%^x) , 



hN,xxx H /, N f /, N H 77-^ {de (hN) Cn,x)x 

{hN)Je{hN) 9s{hNj 



9e 



<AiF{0) = A2, 



and thus 



C^i^' ||-^e,7v||2 < ^2, < t < CT, 



where A2 is independent of A^, e, and ex. Similarly, 



and 



/ dt' {hN) f^e,N,x\\l <As, <t <(T, 

Jo 

[ dt'\\cNJs,N\\l<McN\L<A2kl 0<t<a, 
Jo 



where A2 and A3 are independent of A^, e, and a. By rewriting the evolution equations as 

hN,t = —Je,N,x, ide (^Af) c)j = —Ks,N,x, 

where K^^n = cnJe^u — Qe Q^n) l^e,N,x, we see that there are uniform bounds for Jq dt' || Je,Ar||2 and 
Jq dt' ||i^e,Ar||2, which depend only on the initial data cq (x) and ho (x). We mention the following 
result due Bernis and Friedman |53j . 

[Bernis and Friedman, 1990] Let ipi {x,t) be a sequence of functions that weakly satisfy 
the equation 

^i,t = -Ji,x, Ji = J i^Pi) ■ 

If ft (x, ■) is Holder continuous (exponent \), and if the fluxes Ji satisfy 



dt' W-JiWl < Ai, 0<t<a, 



where Ai is a number independent of the index i and the time a, then there is a constant 
A2, independent of i and a, such that 



lv5i(-,^2) - V5i(-,ti)| < A2\t2 - til 
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for all ti and t2 in (0,cr). 

We now observe that the fluxes Je,N, and K^^j^j satisfy the conditions of this theorem, and thus 

The functions (^t) and Cn {■,t) are Holder continuous (exponent |), for < t < cr. 

We therefore have a uniformly bounded and equicontinuous family of functions {{h^, cn)}n=No+i- 
We also have a recipe for constructing a uniformly bounded and equicontinuous approximate solu- 
tion (/lAT (x, t) , Cat (x, t)), in a small interval (0, a). The recipe can be iterated step-by-step, and we 
obtain a uniformly bounded and equicontinuous family of approximate solutions {{Hn, CN)}N=No+iy 
on (0,To] X Q, for an arbitrary time Tq. Then, using the Arzela-Ascoli theorem, we obtain the 
following convergence result: 

There is a subsequence of the family {{hN,CN)}'^=Ng+i that converges uniformly to a 
limit {h,c), in (0,To] x n. 

In contrast to Ch. [2| this convergence result is uniform and true everywhere in x G i7 (uniform 
strong convergence). We prove several facts about the pair {h,c). 



Let {h, c) be the limit of the family of functions {(/iat, CAr)}^=jVo+i constructed in Steps (5.3.1 )- 
( 5.3.4[ ). Then the following properties hold for this limit: 



1. The functions h (x, t) and c (x, t) are uniformly Holder continuous in space ( expo- 
nent \), and uniformly Holder continuous in time (exponent 

2. The initial condition {h,c) (x,0) = {ho,Co) (x) holds; 

3. {h, c) satisfy the boundary conditions of the original problem (periodic boundary 
conditions ); 

4- The derivatives {h,c)^, {h,c)^, {h,c)^^, {h,c)^^^, and {h,c)^^^^ are continuous in 
the set (0, Tq] x Vl; 

5. The function pair [h, c) satisfy the following weak form of the PDEs, 




dtdx hipt + / / dtdx Js^Px = 0, 
Qtq J J Qtq 




— hx - 

9e{h) g. 

dtdx gs (h) cipt ~^ I I dtdx K^ipx = 0, 




Je = fe (h) hxxx - TTTTY^^ - ^-Z^ [de {h) cl)^. 



Ke = cJe - ge (h) 



c^-c-^(,.(/^)c.). 



where ip (x,t) and ip (x,t) are suitable test functions. 
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The statements (1), (2), and (3) are obvious. Now, any pair {Hn {x,t) ,cn {x,t)) satisfies the 
equation set 




dtdx hN(ft + 




dtdx Je^Nfx = 0, 



9e {h 



N 



9s {h 



/ / dtdx gs (hN) CNijjt + / / dtdx Ki^ j^jiIj^ = Q 



Ke^N = CnJe^N — Qe [h]^) 



- Cn - 



9e (h 



{ge (hN) Cn,x). 



N 



and from the boundedness of the fluxes Js^n and K^^n in L'^ (0,To, (^)), it follows that 

{Je,N, K,^n) ^ ( Je, K,) , weakly in (O, Tq, (fi)) , 

for a subsequence. Using the regularity theory for uniformly parabolic equations and the uniform 
Holder continuity of the (/iat, C7v)'s, it follows that 

The derivatives (/iAr,Civ)j, {hN,CN)^, {hN,CN)^^, {hN.CN)^^^ and {fiNyCN)^^^^ are uni- 
formly convergent in any compact subset of (0,To] x Q. 



Thus, 



Je = fe {h) Kxx - TTTT^^' - "^-^l^y {9e W cl) ^ 



= cJe - Qe (h) 



9e (h) ge 

1 



c — c — 



9eih) 



{9e {h) 



on (0,To] X fi, and therefore Claims (4) and (5) follow. 



5.3.5 Convergence of regularized problem, as £ ^ 



The result in Step (5.3.4) produced a solution (/i^, to the regularized problem. Due to the result 



he (X, t) > /ijnin > B = —kiVL+ \ kfL + 



kfL 



> 0, 



(5.29) 



independent of e, the argument of Step (5.3.4) can be recycled to produce a solution {h,c) to 



the unregularized problem. This solution is constructed as a limit {h,c) = Mnie^o {he, Ce), and 



the results of the theorem in Step (5.3.4) apply again to (/i, c). The result (5.29) applies to h 



constructed as h = limg^o ^e? and thus all the derivatives {h,c)^, {h,c)^, {h,c)^^, (h^c)^^^, and 
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{h, c)^^^^ are continuous on the whole space (0, Tq] x Q and therefore, the weak solution {h, c) is 
in fact a strong one. 



5.3.6 Regularity properties of the solution {h, c) 

Using a bootstrap argument, we show that the solution {h, c) belongs to the regularity classes 



L°° (0, To; /72,2 (Q)) and (0, Tq; H^'^ (fi)). From Step (jsg) it follows immediately that 

Whxh, Wcxh < 

with time-independent bounds. Thus, using Poincare's inequality, it follows that h,c ^ H^'^ (fi) 
and, moreover, 

sup ||/ia;||2, sup ||C2:||2 < OO. 
[O.To] [0,To] 



From Step (5.3.4), it follows that J and Hx belong to the regularity class (0,To;L^ (^)); and 



hence J, Hx G (0,To; iP))- The functions J, /x, and /Za; take the form 

J h hxxx hxh h {jijxC^ ~t~ 'ZhcxCxx^ ? 
= c'^ - c-h~^ {hxCx + hcxx) , 

and 

l^x ^3c 1^ ~t~ h h^Cx h hxx^x ^ ^x^xx ^xxx^ 

respectively. The following results will help us in our demonstration, 

• We shall use the boundedness of h (x, t), 

< /imin < h (X, t) < /imax < OO, 

and the boundedness of c(x, t), ||c||oo < 00. 

• Since fix G (0, Tq; {^)), it follows that fi E L"^ (0, Tq; (f^)), by Poincare's inequality. 

• From this it follows that Cxhx + hCxx is in the class (0, Tq; (fi))- 

• We have the inequality 



||^A''Cx>||2 ^ ^maxll /^ll 00 II '^^ II 2 — ^max||Cx'||2 



^||;U||2 + VZWfixh 
Ij 



and thus hxC^. + hcxCxx £ L"^ (0, Tq; (fi))- 
• Similarly, since J^" dt \\hjjhx\\2 < 00, we have the result h'^Cx + hhxCxx ^ L'^ (0,To; (fi))- 

Now inspection of fi shows that Cxx is in (0, Tq] (fi)), from which follows the result hxxx, Cxxx £ 
(0,To; i^))- By repeating the same argument, we find that Cxx £ L"^ (0,To; (^^))- We also 
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have the result that ||^^Ca;||2 is almost always bounded. To show that hxx G (0,7o;L^ (^))) we 
take the evolution equation for h {x,t), multiply it by h, and integrate, obtaining 



]^ d 



/ dxh^ — 3 / dxh^h^hxx— / dxhxJo= / dxh^h^.^, 
Jn Jn Jo, Jn 



where Jo = hxh ^ + 2h^CxCxx + h'^hxC^. The time integral of the first term on the left-hand side of 
this equation is obviously bounded in time. Let us examine the time integral of the second term, 



/ dt [ dxh^hlhxx < hl^ax i sup \\hx\\l] [ dt\\hxx\\oo 
Jo Jn \[o,To] J Jo 

sup \\hx\\l \ [ dt {L'^\\hxx\\i + \\hxxx\\i) < oo. 

[O.To] J Jo 



— "'max 



The third term on the left-hand side is despatched with in a similar way, so that Jq^° c?t||/ia;x||2 < co- 
We have now shown that 

Using this result, together with the previous facts gathered together in this section, it is readily 
shown that 

h^ceL" {0,To;H^'' (Q)) , 
And finally, using this result, it follows that 

For example. 



/ dt [ dxhlcl^ < (sup\\hx\\l] [ 
Jo Jn \[o,To] J Jo 



dt 1 1 Cxx 1 1 oo 



< ( sup ||/la;||2 ) / dt [L ^\\Cxx\\l + \\Cxxx\\lY < OO. 
\[0,To] J Jo 

This bound, together with the result h^hx G (0,To; (^)), implies that 

hlcxeL^ {0,To;L^ (Q)) . 

Similarly, 

[ dt j dxhlcl^<\ sup\\hx\\l \ [ ci^ (l/"^||Ca,a,||i ||Ca;^^||l)^, 
Jo Jn \[o,To] J Jo 

Hence hxCxx G (0, Tq; (^)), and by symmetry, Cxhxx G (0, Tq; (f^)). Since 
Hx — (Sc 1^ Cx ~\~ h h^Cx h hxx^x h hxCxx ^xxxj 
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is in (0,To; (f^)), it follows that Cxxx is in this class too. A similar argument holds for hxxx- 
Thus, the solution {h, c) belongs to the following regularity class: 

{h, c) G (0, To; H^'^ (Q)) n (O, Tq; H^'^ (Q)) n C^-i {Q x [0, Tq]) . (5.30) 

Extra regularity is obtained by writing the equation pair as 

dh 

— + h^hxxxx = -^h^hxhxxx + '^i + '^2 = V (a;, t) , 
dc 2 , , , / \ 

^ + Cxxxx = -J^hxCxxx + Ipl + ^2 = W [X, t) , 

where 

/ c^^llv^ib < sup||c^^||2 / dt\ui\, ui e {[0,t]) , 

Jo \[0,r] J Jo 

/ dt\\ilJi\\2 < \ snp Wcxxh ] / dt\u2\, U2 e {[0,t]) , 

Jo \[0,r] J Jo 

for any r G (0, Tq], and where ip2 and 4'2 belong to the class (0, Tq; (f^))- By multiplying the 
height and concentration equations by hxxxx and Cxxxx respectively, and by integrating over space 
and time, it is readily shown that 

(/i,c)GL- {0,To,H''' (fi)), 

and hence 

from which follows the regularity result 

(h, c) G (0, To; H^'^ (^])) n (O, Tq; H^'^ (fi)) n Ci'S (^] x [0, Tq]) . (5.31) 



5.3.7 Uniqueness of solutions 



Let us consider two solution pairs {h, c) and {h', c') and form the difference {6h, 6c) = {h — h',c — c'). 
Given the initial conditions {Sc (x, 0) , Sh (x, 0)) = (0, 0), we show that {6h, 6c) = (0, 0) for all time, 
that is, that the solution we have constructed is unique. We observe that that the equation for 
the difference 6c can be written in the form 



d 



(5.32) 



where 6(f (x, t) G (0, Tq, (^)), and where 6ip {6c = 0) = 0. Thus, using the semigroup theory 
of Sec. 2.5, we find that Eq. ( 5.32[ ) has a unique solution. Since 6c = satisfies Eq. (5.32), and 
since 6c (x, 0) = 0, it follows that 6c = for all times t G [0, Tq]. 
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It is now possible to formulate an equation for the difference Sh by subtracting the evolution 
equations of h and h' from one another, mindful that Sc — 0. We multiply the resulting equation 
by Shxx and integrate over space, obtaining the result 



-A [ dx6hl+ [ dxhHhl^= [ dx6hlJhx{h-' + h''cl) 
"J Jn Jo, Jn 

- dx {h^ - h'^) {h'^^^ + 4:CxCxx) 



Jn 



Using the lower bound on h {x, t) > /imin > and Young's first inequality, this equation is trans- 
formed into an inequality, 

1 ^ ^ + dx 6hl^^ - '^^ ^'^ '^^^^^ ^ Z ^^"^ ^ ^'^^^^ ^ 

+ K, [ dx5hl^ + ^ [ dx {h'-h'')\h',^, + ACxCxxf 

Jn 4/^2 Jn 

+ kJ dx 5hl,, + ^ [ dx [K {h-' - h'-') + hxcl {h' - h")] ' , 
Jn ^1^3 Jn 

where ki, K2, and K3 are arbitrary positive constants. By choosing ki + K2 + — h^^^, the 
inequality simplifies, 

l4- [ dxShl<-^ [ dx6hl(h-' + h'^cl? 



+ [ dx {h^- h'^Y {h'^^^ + ACxCxxf 

4/^2 Jn 

^ [ dx [K {h-' - h'-') + hxcl {h' - h'^)] ' , 
\^i^3 Jn 



+ , 



which in turn reduces to 

-1 , L/2 2\2 



y dx5hl< J dxShl {h-^ + h'^lf 

) +{K + hxciy 

Jn L 

where k, is another positive constant. We integrate over the time interval [0,7], and use the fact 
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that ||(5^a;||2 (0) = to obtain the string of inequahties 



2k sup ||5/i2:||2 (t 

Te[o,T] 



< [ dt [ dx6hl {h-^ + h'^ciy 
Jo Jn 



< / dt\\SK\\l\\h-^ + h'^c 



1 , U'2„2u2 
a; II oo 



+ / dt dx6h 



dx 



< sup \\6K\\1{t) I dt\\h-^ + h''^c: 

re[0,T] 



+ / dt\\6h\\ 



1 , U'2„2\\2 

a; II oo 



iKxx + ^C^C^xf + {K + hxCl)' 



+ sup \\5h\\l^{T) I dt I dx [Kxx + Ac^c^xf + [h'^ + Kclf 

re[0,T] Jo JQ. ^ 

The Poincare inequahty can be combined with the one-dimensional differential inequalities dis- 
cussed in Ch. |2]to yield the relation ||/||oo < ^p||/a;||2, where / is some mean- zero function and Kp 
is an /-independent constant. We therefore arrive at the inequality 



2k sup \\6K\\1{t)< sup \\5K\\1{t) [ dt\\h'' + h'^cl\\l 

re[0,T] Te[0,T] Jo 

+ Kp sup \\Shr,\\l{T) / dt [\\h'^^^ + 4cxCxx\\l + \\h'^ + h^cllll) . 
tg[o,t] Jo 



Using the results of Sec. 5.3.6, it is readily shown that + /I'^c^ G L"^ (0,T; {^)), and that 
the functions h'^^^ + iCxCxx and h'^ + hxcl belong to the class (0,T; (f^))- By choosing T 
sufficiently small, it is possible to impose the inequality 



1 

2k 



I dt\\h~' + h"cl\\l + Kl [ dt {\\h',^^ + AcxCxxWl + \\K + KclWl) 
Jo Jo 

which in turn forces sup^g^ y] ||5/ik||2 = 0, and hence the solution is unique. 



< 1, 



5.4 Summary 

In this chapter we have switched focus from the passive to the active tracer. The jump in complexity 
is significant, and involves the Navier-Stokes Cahn-Hilliard equations. We focus on the case where 
the active system forms a thin layer on a substrate, and obtain simplified thin-film Stokes Cahn- 
Hilliard equations that are valid when the system's vertical gradients are small compared to the 
those in the lateral directions. We analyze these equations in one lateral dimension, using the tools 
developed in Ch.[2| and prove existence, regularity, and uniqueness results for the equations. While 
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this analysis yields many precise results, it is unable to provide qualitative information about the 
solutions, in particular the effect of the backreaction on the concentration and free surface height, 
and we therefore turn to numerical simulations in the next chapter. 



Chapter 6 



Nonlinear dynamics of phase separation 
in thin films: Qualitative features 



6.1 Overview 



We take the thin-film Stokes Cahn-HiUiard (SCH) equations in two and three dimensions as a 
starting point and investigate their quahtative features. By inspection of Eqs. (5.16), we can 



write down a three-dimensional set of SCH equations for the free surface h(x,t), and for the 
phase-separating concentration field c{x,t), and it is these we shall study in this chapter. 



dh 
di 



+ Vx- J = 0, 



Ft 



[ch) + V± ■ (cJ) = V± ■ (/iV±/x) 



where 



J = U2v,r 



1 

c 



h 



{hV±c) . 



(6.1a) 
(6.1b) 

(6.1c) 
(6.1d) 



We have used the derivative operator in the lateral directions, Vj 
the following nondimensional constants. 



ho 



C 



upD 



2' 



{dx, dy), and have introduced 
ho 



where e is the ratio of the horizontal to the vertical lengthscales, assumed to be small. In Sec. 6.2 



we 



shall study the equihbrium version of these equations and, inspired by Sec. 3.6, we shall reduce this 



problem to a boundary-value problem. In Sec. 6.3 we shall carry out simulations of the equations 
in three dimensions and study mechanisms for controlling the phase separation. Crucially, we find 
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that these mechanisms depend on the value of the backreaction strength, and we discuss this fact 
in the context of microfabrication apphcations. 



6.2 Equilibrium solutions of the model equations 



In this section we discuss equihbrium solutions of the model equations (6.1 ). We use this knowledge 



to show that in the presence of a Van der Waals regularizing potential, the system exhibits the 
spinodal instability with an ultimate state of phase separation and a free-surface height slaved to 
the concentration field. For simplicity, in this section we shall restrict to one lateral dimension. 

Let us set = A/h^ and F = constant. We first of all investigate the stability criteria for the 



constant solution (/io,Co) to Eq. (6.1) by identifying small deviations away from the base state, 
6h oc e^ht^ikx^ ^ ^cTct^tkx^ Insertion of this ansatz into Eq. (6.1) gives 



hie 



C ' 
1) e 



A 

K 



(6.2) 



2l4 



cik 



The system is stable to perturbations in the height field if A < 0. For A > 0, the height field 
will rupture in finite time [301 llU2j . For |co| < 1/v^ however, the spinodal instability is accessible 
independently of the sign of A, which suggests that a critical mixed state will phase separate in a 



manner similar to the classical Cahn-Hilliard fluid, as described in Sec. 2.2 



We are interested in the interplay between concentration gradients and the height profile at late 
times and so we set A = — |y4| and study the system in equilibrium, in one dimension. By inspection 
of Eq. (6.1), the equilibrium conditions are ^ = constant, u = which gives equations 



1 d^h 



C'\A\ 1- 



1 



(6.3a) 



d'^c 3 Idh dc 

dx"^ hdx dx' 



(6.3b) 



where we have enforced the boundary conditions h{±oo) = 1, /i (±oo) = and have rescaled 
lengths by C^. 



For the case C = oo and p = r/C'^\A\ <^ 1, Eq. (6.3) has an asymptotic solution. We find the 
/i-equation 



Hence, the c-equation is 



h= ll + p 



1 /'^2 
4 



-1/3 



(6.4) 



(C3 - C) 



(6.5) 
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For small p, the solution is c = tanh (x/ v^) + O (p) and hence 



h = 1 — ^psech^ 



/ X 

V71 



+ 0(p^) 



p < 1. 



(6.6) 



Thus, in this limiting case, the height profile is approximately constant {h = 1) except in the 
transition region of the concentration field, where the height profile possesses a valley. 

For nonzero surface tension, two parameters characterize the problem, for now the equations are 



CCi \A\ I 1 



1 



Cr 



dx"^ 



1 dh dc 



c — c — 



hdx dx' 



(6.7a) 



(6.7b) 



the two independent parameters being CC"^ \A\ and Cr. In Fig. 6.1 we present numerical solutions 
exhibiting the dependence of the solutions on these parameters. As before, the height field possesses 
peaks and valleys, where the valleys occur in the transition region of concentration. While the 
valley increases in depth for large r, rupture never takes place. This result follows from the 
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Figure 6.1: Equilibrium solutions of the thin-film equations in one dimension, obtained from solving 
a boundary-value problem with parameter values C = C'^\A\ = 1 and r = 0.1, 1, 10, 50. In (a) the 
valley deepens with increasing r although the film never ruptures, while in (b) the front steepens 
with increasing r. In (c) we plot the forces F^^^ and -FVdw for C = = |A| = 1 and r = 50 to 
verify that they have opposite sign. 



inequality h" (0) > 0, since a; = is a local minimum. Thus, from Eq. (6.7) we obtain the bound 

'\ {c{of-iY + lc'{oy 



< [h 



(6.8) 



The repulsive Van der Waals potential therefore has a regularizing effect on the solutions. We 
carried out a further numerical integration to find the general dependence of the magnitude of the 



dip on the parameters of Eq. (6.7). There is no clear scaling law for this dependence, although 



the approximate rule hj^^ ~ r/C^|y4| is evidenced by Fig. 6.2, a relationship that is exact when 
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1 1 I 600 




(c) (d) 

Figure 6.2: The dependence of dip magnitude /imin on the problem parameters. Subfigures (a) 
and (b) show the dip magnitude's dependence on the parameter C^l^l; Subfigures (c) and (d) give 
the behaviour of the dip magnitude as a function of r. In both cases, the approximate relation 
^mfn ~ ''^/^nl^l discemablc. 



C = oo, as seen in Eq. (6.4). Thus, we have characterized the equilibrium solution as fully as 



possible. We now consider the physical meaning of these results. 

The formation of the valley in the height field has the physical interpretation of a balance between 



the Van der Waals and backreaction effects. From Fig. 6.1 (c) we see that the backreaction or 
capillary force Fcap = —rh~^dx \h{dxcf~\ and the Van der Waals force -FVdw = \M'^xh~'^ have 
opposite sign. The repulsive Van der Waals force acts like a nonlinear diffusion and inhibits 
rupture, and therefore Fcap promotes rupture, a result seen in experiments [39]. The valley in the 
height field therefore represents a balance between the smoothening and rupture-inducing effects. 



We compare the no-rupture condition in Eq. (6.8) with the analytical results obtained in Ch. 5 



In terms of the physical parameters of the system, the no-rupture condition of Ch. [5] is 




/i^in > B = J2CLiFo + Fi\A\) \ . 1 > 0, 

where Fi = \ j^^ dx [h (x, 0)]~^ ^ 0, and Fq = F (0) — Fi. The function 5 C) has no exphcit r- 
dependence: although Fq depends on r, it is possible to find initial data to remove this dependence. 



6.2. Equilibrium solutions of the model equations 



118 



We show a representative plot of B {\A\,C) in Fig. 6.3 Although a comparison between Fig. 6.2 



and Fig. 6^ is not exact, since the boundary conditions and domains are different in both cases, 
we see that the shape of the bound in Fig. |6.2 is different from that in Fig. 6^ Since the bound in 



Fig. 6.2 is obtained from numerical simulations, and is intuitively correct, we conclude that it has 



the correct shape and that the bound of Fig. 6.3, while mathematically indispensable, is not sharp 
enough to be useful in determining the parametric dependence of the dip in free-surface height. 



0.02 




Figure 6.3: A typical plot of B {\A\, C) for Fq = Fi = ^ and C = 1. This theoretical lower bound 



has a different shape from those in Fig. 6.2, which suggests that while B {\A\, C) plays an important 



role in the analysis of the model equations, it does not capture the physics of film thinning. 



The procedure of identifying a phase-separating instability and studying non-trivial equilibrium 
solutions is recurrent in Cahn-Hilliard dynamics, and is discussed in Ch. |2} In that case, the non- 
trivial one-dimensional equilibrium solution hints at the late-time solution in several dimensions. 
Indeed, the one-dimensional antikink solution suggests that the late-time solution in several dimen- 
sions consists of domains with transition regions of the same width as the antikink width. Thus, 



we expect the equilibrium solution in Fig. 6.1 to reflect the late-time structure of a general phase- 



separating fluid in a thin film. This is found to be the case in one-dimensional, time-dependent 



numerical simulations, of which Fig. 6A_ is an example; here the concentration forms domains 
c ~ ±1, separated by smooth transition regions, while the free-surface height forms peaks and 
valleys, with valleys occurring in the transition regions. The flatness of the peaks depends on the 
Van der Waals parameter |y4|, with large 1^41 producing flat peaks, while the depth of the valleys 
depends on the value of the backreaction strength r. The periodicity in Fig. 6.4| is due to the 



periodicity of the initial conditions, in particular the period-four initial condition for the height 
field. Purely random initial conditions give rise to a similar but aperiodic pattern of domains, 
peaks, and valleys. 

We have found by one-dimensional calculations, that the ultimate state of the equation system 



Eq. (1.1) is a concentration field comprising domains, with a height profile that dips at domain 
boundaries. In order to pursue this claim further, and to examine the dynamics of domain growth, 
we turn to simulations in two dimensions. 
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Figure 6.4: (a) A numerical solution of the SCH equations in one dimension, with C = —A = 1 
and r = 1 (large backreaction); (b) The evolution of the free-energy functional. The free energy 
decays monotonically, in agreement with the theory, except at late time, where there is a slight 
loss of monotonicity due to numerical error. 



6.3 Solutions in two dimensions 



In this section we outline a numerical method for the two-dimensional thin-film Cahn-Hilliard 
equations and apply it in two distinct cases: constant surface tension only, and spatially-varying 
surface tension (or applied tangential stress). In each case, the regularizing Van der Waals potential 
is included to avoid film rupture. The introduction of appropriate structure functions enables us 
to measure the lengthscales of the concentration morphology. 

To measure the typical scale of the domains we introduce a structure function s {k, t) in what 
follows. Because the concentration field is dependent only on the lateral coordinates x and y, we 



work in two dimensions only. Thus, as in Sec. 3.2 a measure of the correlation of concentration 



between two neighbouring points in Fourier space is given by the expression 



cix + x') c (x) — c ^ 



where L is the box size, x = {x,y) and where (■) denotes the spatial average. We normalize this 
function and compute its spherical average to obtain the structure function 

where the spherical average (p (k) any function (p (k) is defined as 



<p{k) = —i de<p{k). 







Hence is the spherical average of the Fourier coefficient c^. Thus we identify the dominant 
lengthscale R (t) with the reciprocal of the most important fc- value, defined as the first moment of 
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the distribution s{k,t), 

In performing the spherical average, we have assumed that the system is isotropic. If the isotropy 
is broken, for example, by the presence of an external shear flow, then a vector-valued measure of 
the domain lengthscales is given by 

Jd^kkSik,t) 

Using the initial conditions and the scale analysis we have outlined, we carry out the simulations 
for a variety of external conditions. 



Numerical method 



To simulate the thin-film SCH equations (6.1), we use a semi-implicit finite-difference algorithm 



in two dimensions, although here we sketch the method in one dimension only. The height field is 



updated first, using semi-implicit finite differences |103] . The height field of Eq. (1.1) is discretized 
in time as follows, 

{h-^'Y^4z^ + suh'',cn, (6.10) 



At 



3C 



where 



Sh {h, c) = K- 



.dh d ( 1 d^h _ \A\ 
dx dx \ C dx^ 



3 dx'^ 



=po 



1 



+ 



r d 



d 



3 (9x 1 dx 



h(—\ 
\dx J 



contains third-order derivatives and lower. Equation (6.10) is solved for /i""*"^ using the Newton 



method for root-finding. Using the updated value of the height, we update the concentration, 
making use of the following discretization in time. 



At 



n4 n+l 



(6.11) 



where 



Sc{h,c) 



dc 1 dh da d'^ 

-u 1 — H 

dx h dx dx dx"^ 



1 dh dc 
h dx dx 



contains third-order derivatives and lower. Here u is the velocity and /i is the chemical potential. 



defined in Sec. 5.2 Equation (6.11) is solved using a spectral method |5]. 



We solve the equations on a 128 x 128 square grid with a timestep of At = 5 x 10~^, and the 
results do not change markedly upon increasing the resolution or decreasing the timestep. Following 
standard practice [161 [20l [71] , the Cahn number is taken to be a/ Ax'^/At so that the transition 
region is just barely resolved. The boundary conditions are periodic in the lateral directions. 
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The initial condition for the height field is a sinusoidal perturbation while the concentration field 
is given a random perturtion about zero. In contrast to the usual Cahn-Hilliard equation [71] 
however, the zero spatial mean of the concentration field does not persist, since it is the quantity 
J dPxc (a?, t) h (a;, t), not / (Pxc (a?, t) that is conserved. 



Constant surface tension 



We let r = constant, and study the velocity field 



u 



Vj 



\A\ 



h 



(6.12) 



We carry out simulations according to the scheme discussed above and find that the system evolves 
towards towards its late-time configuration in the following way. At early times, domains of like 
concentration start to form, while the initial configuration of the free surface is destroyed due 



10 



10 




Figure 6.5: Log-log plot indicating the time dependence of the typical domain size R{t). 
dependence is close to the Lifshitz-Slyozov growth law, R ~ t^/^. 



This 



to nonlinear diffusive effects in the height equation. At later times, the domains continue to 
form and then to grow, while the free-surface height becomes slaved to the concentration field. 
That is, the spatial and temporal variations of the free surface depend on the concentration level. 
The free-surface height possesses local maxima in regions where domains of concentration have 
formed (peaks) while at the transition between domains of different binary fiuid components, the 



free-surface height decreases markedly (valleys), as in the one-dimensional simulations of Sec. 6.2 



The location of the surface peaks over the concentration domains has been seen in experiments 
involving polymeric thin films jlD4] . As in the case of ordinary Cahn-Hilliard dynamics, the late- 
time configuration selected is the one that tends to minimize the free energy. Now the free energy 
has the form 



F[cM 



1 

2C 



Vh\' + 



(fx + 



C2 



h 



[o,Lr 



Vcl 



with -F < 0, and the dynamics outlined do indeed decrease this free energy. 
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We characterize the growth of the domains by examining the function R (t) in Fig. 6.5 After tran- 
sient effects have passed, the lengthscale R{t) grows approximately as t^^^, which is the Lifshitz- 



(a) 




(b) 



Figure 6.6: (a) Concentration field at t = 5,000, r = 0.01, C = 1, and A = —10; (b) The height 
field for the same parameter values. Note that the height field is slaved to the concentration field. 



Slyozov growth rate for phase separation in the absence of flow |1H]. Several numerical simulations 
with different values of \A\, C, and r, indicate that the scaling law is independent of these parame- 
ters. A better estimate of the growth law could be obtained by increasing the resolution, although 
this is computationally costly. However, the qualitative features of the system are well described 
by low-resolution simulations. The recovery of the growth law R (t) ~ t^^^ demonstrates that 
the Lifshitz-Slyozov mechanism of evaporation-condensation dominates over other effects, such as 
viscosity-driven phase separation [1]. This is not surprising, since the timescale that determined 
the asymptotic ordering in the derivation of Eq. (1.1) was the diffusion time tdis ~ L'^/D. In this 



sense, our model equations demonstrate the importance of variations in the free-surface height, 
rather than hydrodynamics. At very late times, the scaling rules are spoiled by finite-size effects. 



Variable surface tension 



We let r = Fosinfcx, Fq is a dimensionless amplitude, k = (27r/L)n = kon is the spatial scale 
of the surface tension variation, and n is an integer. Then the velocity that drives the system 
becomes 



u = \h (/cFo cos kx, 0) + \h?' 



Vj 



(6.13) 



By inspection of Eq. (5.3b), this velocity field may also be obtained by imposing a shear stress r 
at the surface, provided r = V_lF. The imposed shear stress acts as a driving force. 

This choice of velocity field leads to control of phase separation in the following manner. For small 
values of the backreaction strength, with r — > 0, the height field quickly aligns with the surface 



tension profile as in Fig. 6.9, since the strong effect of the Van der Waals diffusion destroys the 
unforced part of h{x,t). At the same time, the concentration field begins to form domains. At 
later times, when kx{t), ky [t) ~ k, the domains align with the gradient of the forcing term. The 
growth of the domains continues in this direction and is arrested (or slowed down considerably) 
in the direction perpendicular to the forcing. The domains are string-like, with kinks occurring 



6.3. Solutions in two dimensions 



123 




(a) (b) (c) (d) 

Figure 6.7: Domain structure at t = 1000, t = 3750, t = 7500, and t = 15000 for C = 1, A = -10, 
and r = 0. The surface tension gradient is parallel to the arrow and F = Fgsin (kx), Fq = 20 and 
k = Ako- For r = 0, the domains align along the arrow. 



along lines where F {x,y) is minimized, as evidenced by Fig. 6.7, The decay of fc^,. and ky is shown 



in Fig. 6.10 It is not clear whether the decay of {k^, ky) is arrested or continues slowly, and we do 



not report the decay rate here. 

For moderate values of the backreaction strength with r ~ 0(1), the height field again assumes a 
profile aligned with the surface tension, while domains of concentration now align in a direction 
perpendicular to the forcing gradient. Domain growth continues in the perpendicular direction and 

HI mmm 



(a) 



(b) 



(c) 



(d) 



Figure 6.8: Domain structure at t = 1000, t = 3750, t = 7500, and t = 15000 for C = 1, A = -10, 
and r = ^. The surface tension gradient is parallel to the arrow and F = Fq sin (kx), Fq = 20 and 
k = 4ko. The system reaches a steady state where the domains align perpendicular to the arrow. 



is arrested in the direction of the driving-force gradient. A pattern of string- like domains emerges, 
with domain boundaries forming along lines where both a (x, y) and h {x, y, t) are maximized. 
Eventually, the domain boundaries align perfectly with the surface tension maxima, as evidenced 
in Fig. |6.8 



The control of phase separation therefore depends crucially on the backreaction. This result is 
amplified by the existence of a no-rupture condition only for the r = case, where there is no 
backreaction. This condition relies on the alignment of the height and surface tension profiles, 
which is exact only when the backreaction is zero. Then at late times, the system evolves towards 
equilibrium and is described by the steady state V ■ [^h'^V±T + |/i^V_L {C^^V^h + |y4|/i"^)] = 0, 
which by the alignment property reduces to the one-dimensional equation 



1 

h -h — 

dx ^ dx 



1 d^h 
C dx' 



+ 



\A\ 



constant. 



(6.14) 
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Figure 6.9: The free-surface height for r = and t = 15000 ahgns with the apphed surface tension. 
The height field at t = 15000 for r = | is similar. 

By multiplying both sides of the expression by h, differentiating and then evaluating the result at 
Xq, a coincident minimum of height and surface tension, we obtain the condition 



[hixo)\ 



\A\h"ixo). 



Since xq is a minimum of height, h" {xq) > 0, which prevents h (xq) from being zero. On the other 
hand, for r and Fq sufficiently large, the alignment of height and surface tension profiles is not 
exact, the reduction of the steady-state condition V ■ [^h'^'V±T + lh^'V± [C'^V^h + {Alh'^)] = 



to Eq. (6.14) is no longer possible, and our calculation fails. In that case, simulations show that 



the film ruptures in finite time. 

We have outlined, by numerical simulations and calculations, three possible outcomes for the phase 
separation, depending on the backreaction strength r. For r -C 1, the concentration forms string- 
like domains, aligned with the applied force. For r = O (1) the concentration forms domains 
that align perfectly in a direction perpendicular to the applied force, while for r ^ 1, the forcing 
causes the film to rupture. Loosely speaking, the different regimes occur because for r = the 
concentration field couples to the gradient V_lF, while for r ~ O (1), the concentration field couples 
to the patterned surface tension F. Thus, the interfacial tension or backreaction can therefore be 
tuned to achieve the desired outcome. 



6.4 Summary 

We have derived a thin-film model of phase separation based on the NSCH equations, in which 
the reaction of concentration gradients on the fiow is important. We have used this model to 
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(a) (b) 

Figure 6.10: Time dependence of the dominant spatial scales and ky for (a) r = 0, where the 
secular behaviour of k^ and ky has a small drift whose rate we do not report; (b) r = |, where 
kr ^ 2 and k,, 0. 



give a qualitative picture of the features of phase separation in real thin films, in particular the 
tendency of concentration gradients to promote rupture in the film, and to produce peaks and 
valleys in the free surface that mirror the underlying domain morphology. We have found that 
in the presence of a unidirectional sinusoidal variation in surface tension, the strength of the 
backreaction determines the direction in which the domains align. This result could prove useful 
in microfabrication applications where control of phase separation is required [37] . 

Because the lubrication model suppresses vertical variations in the concentration field, we are 
limited to the case where the binary fluid components interact identically with the boundaries at 
the substrate and free surface. However, the model quite generally gives an accurate description 
of surface roughening arising from Van der Waals forces. More detailed models based on this 
approach, involving different boundary conditions that better reflect wetting behaviour and a 
concentration-dependent Hamakar coefficient, could capture a wider range of thin- film behaviour, 
although the tradeoff could be a loss of mathematically interesting features, such as the Lyapunov 
functional. 



Chapter 7 



Singular solutions of aggregation 
equations 



7.1 Overview 

In this chapter we introduce a generahzed gradient flow equation that serves as a bridge from 
the Cahn-Hilhard theory we have discussed, to models of aggregation that are currently used 
in nanoscale physics. We then focus on a particular aggregation model that describes a nonlocal 
magnetic system and admits singular solutions that emerge from smooth initial data. The magnetic 
system comprises a Gilbert-type equation for an orientation vector, and a scalar density equation. 
We study the uncoupled Gilbert-type equation and find that the evolution depends strongly on 
the lengthscales of the nonlocal effects. We then pass to a coupled density-magnetization model 
and perform a linear stability analysis, noting the effect of the lengthscales of nonlocality on 
the system's stability properties. We carry out numerical simulations of the coupled system and 
find that singular solutions emerge from smooth initial data. The singular solutions represent a 
collection of interacting particles (clumpons). By restricting ourselves to the two-clumpon case, we 
are reduced to a two-dimensional dynamical system that is readily analyzed, and thus we classify 
the different clumpon interactions possible. 



7.2 Gradient flow, variable mobility, and nonlocality 

In this section we introduce a generalized aggregation equation that includes the Cahn-Hilliard, 
Keller-Segel, and (scalar) Holm-Putkaradze equations as special cases. 

Recall that the Cahn-Hilliard equation for an order parameter ■0 is derived from the gradient of a 
free-energy functional in n dimensions. 




126 



7.2. Gradient How, variable mobility, and nonlocality 



127 



where /tsm (V') is a potential with two stable minima. In the previous chapters, we chose the 



double-well potential /tsm (tp) = \ (V^^ — 1) , and this choice simplified the analysis of the model. 
However, any potential with two stable minima and one unstable maximum will do. Therefore, in 
this chapter, we choose to work with the so-called double-obstacle potential 

/tsmW =eW^W+V'(l-V'), =V'l0gV' + (l-^)l0g(l-V'), 

where ^ is a positive parameter and the function W is minus the system's entropy density. Recall 
that the quartic double well was invariant under the component relabelling transformation ip — > 
The logarithmic double well function has a similar symmetry: it is invariant under the 
relabelling transformation —>■ 1 — ip. We can now write down the Cahn-Hilliard equation in 
nonlocal form, 



CH 



e f d^xW{,lj)+ [ d^xi:{x,t) [ ryKcu{\x-y\){l-ijiy,t)) 
Jq Jn Jq 



where Ken {\x — y\) is the kernel whose Fourier transform is 



K(k) = l- ^^k^ 



It is possible to prescribe a different kernel: other choices can lead to anisotropic or nonlocal interac- 
tions between fluid particles. Let us therefore examine the possibility of nonlocal interactions, and 
leave the kernel unspecifled. We shall, however, restrict to central interactions, K [x) = K{\x\). 
The free energy functional is thus 

F[ij]=i f <rxW{^)+ [ <rx^{x,t) [ ryK{\x-y\){l-^P{y,t)). (7.1) 
Jci Jvl Jn 

As discussed in Ch. |2| the simplest rotationally-invariant, mass-preserving evolution equation that 
decreases this free energy is the gradient flow equation 

where M (tp) is a nonnegative function called the mobility and fi = 6F/6ip is the chemical potential. 
This evolution decreases the free energy: 



dF 



d'^xM (^) |v/i|^ 



The mobility function controls the strength of the particle diffusion and depends on the local 
density of the fluid. If the evolution is to respect the component relabelling symmetry, then the 
mobility must be a function of the combination ip [1 — ip). The simplest possible choice of mobility 
is thus 

M {ip) = {I - ip) . 

The evolution of the system preserves the positivity of this function because the energy cost of 



attaining the either of the states ■?/' = 0, 1 is inflnite, as expressed in Eq. (7.1). 
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The nonlocal Calin-Hilliard concentration for the concentration ip {x,t) is thus 



~dt 



V- 



^ (1 - V^) V (^^W (i/j) - ^ d-yK {\x-y\)^ {y, t) 



(7.2) 



By specifying a particular kernel, or by modifying the mobility, we recover a number of models of 
aggregation. In particular. 



• If the Fourier transform of the kernel K \s 1 — ^7^^, we recover the Cahn-Hilliard equation 
with the double-obstacle potential and variable mobility; 

• If the Fourier transform of the kernel K is (A^A;^ + 1) ^, where A is the lengthscale of the 
potential energy d"'yK {\x — y\) {1 — ip {y, t)), that is, if K is the kernel of the Helmholtz 
operator 1 — A^A, we recover the regularized Keller-Segel equatioiQ 

• If i^' is any nonlocal kernel, ^ = 0, and the mobility {1 — tp) is replaced with ip [1 — H *ip), 
where H * ijj is the smoothened density 

{H*ij){x,t)= [ d''yH{\x-y\)ij{y,t) 
Jn 

we obtain the scalar Holm-Putkaradze model 1381. 



The scalar Holm-Putkaradze model was introduced in and describes the evolution of a 

single scalar density. In particular, given the (negative) free energy 



"HP 



I / d"xip{x,t){l-a^A) ^'4j{x,t) 
Jn 



the model admits singular solutions that emerge from smooth initial data. The singular solutions 
satisfy a set of ordinary differential equations. The model can therefore be used to describe self- 
assembly in nanoscale physics, in which atoms or molecules aggregate to from larger, mesoscale 
structures with particle-like properties. These particles are governed by the ordinary differential 



equations to which the aggregation equation reduces. The aggregation equation (7.2) that we have 
formulated thus provides a bridge from the Cahn-Hilliard theory we have discussed, to this more 
recent application in nanoscale physics. The authors Holm, Putkaradze, and Tronci have extended 
this model to a coupled scalar- vector system in [62] , although they simply formulate the equations, 
and do not examine the properties of the model in any detail. Therefore, using the numerical tools 
that we have developed, and our knowledge of aggregation through energy minimization, we turn 
to this coupled scalar-vector system, and provide qualitative and quantitative descriptions of its 
behaviour. 

We treat the initial state of the system as a continuum, a good approximation in nanophysics 
applications |105j . One realization of this problem is in nanoscale magnetic systems, in which 
particles with a definite magnetic moment collapse and form mesoscale structures, that in turn 



*The regularized Keller-Segel equation is discussed in Appendix B 
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have a definite magnetic moment. Thus, in this chapter we refer to the orientation vector in our 
continuum picture as the magnetization. 



7.3 The nonlocal Gilbert equation 

In this section we study a magnetization equation that in form is similar to the Gilbert equation, 
that is, the Landau-Lifshitz-Gilbert equation in the over-damped hmit \106\ 1107] . The equation 
we focus on incorporates nonlocal effects, and was introduced in [62] . We study the evolution and 
energetics of this equation, and examine the importance of the problem lengthscales in determining 
the evolution. 

We study the following nonlocal Gilbert (NG) equation 

dm f 6E\ , 

m X Ui„ X — , (7.3) 



dt \ 5m ^ 

where the magnetization vector m, its smoothened version /x^, and the variational derivative 
6E/6m are functions of the form 

m, fim, SE/6m : (QcR) x (M+ U {0}) R\ 

The smoothened magnetization is defined as 

(1-/3292) -'m, 

and 6E/6m is prescribed as follows, 

(1 - a^dl) m. 



Sm 

The smoothened magnetization and the force 6E/6m can be computed using the theory of 
Green's functions. In particular, 

{x, t) = dyHp {x-y)m {y, t) := Hp*m {x, t) . 
Jn 

Here * denotes the convolution of functions, and the kernel Hp (x) satisfies the equation 

l-f3'^)H,{x) = S{x). (7.4) 



The function S (x) is the Dirac delta function. Equation (7.4) is solved subject to conditions 
imposed on the boundary of the domain Q. In this chapter we shall work with a periodic domain 
Q = [— L/2, L/2] OT fl = [0, L], although other boundary conditions are possible. Since 6=*=**^^, k G 



is a simultaneous eigenfunction of the operators (1 — a^S^) ^ and (1 — P'^d^) \ equation (7.3) has 
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a family of nontrivial equilibrium states given by 

f^eq (x) = mo sin {kx + 0o) , 
where mo is a constant vector, k is some wavenumber, and 0o is a constant phase. The derivation 



of this solution is subject to the boundary conditions discussed in Sec. 7.4 



By setting /5 = and replacing (1 — ce'^d'^) ^ with —d"^, we recover the more familiar Landau- 
Lifshitz-Gilbert equation, in the overdamped limit |106j . 



dm 



—m X m X 



dx"^ 



(7.5) 



Equation (|7.3|) possesses several features that will be useful in understanding the numerical simu- 

(7.6) 



lations. There is an energy functional 

E{t) = ^ / dxm ■ (1 - a'^diy^ m, 
Jn 

which evolves in time according to the relation 



dE 
~dt 



dx 



Hm- [l- a'^dl) ^ m m ■ (l - a^d^) ^ m 

dx {nm ■ n^) {l — a'^dl) ^ m 



dx 



.2q2\-1 



(7.7) 



This is not necessarily a nonincreasing function of time, although setting (3 = gives 



(5=0 



dx 



m ■ (l — a'^dl) ^ m 



dxm 



dxm' 



(l - a^dl) ^ m ' (cos^ <^ - l) < 0, 



(7. 



where cp is the angle between m and (1 — a^d^) ^ m. In the special case when /? — 0, we therefore 
expect E (t) to be a nonincreasing function of time. On the other hand, inspection of Eq. (7.7) 
shows that as a — 0, the energy tends to a constant. Additionally, the magnitude of the vector m is 
conserved. This can be shown by multiplying Eq. (7.3) by m, and by exploiting the antisymmetry 
of the cross product. Thus, we are interested only in the orientation of the vector m; this can 
be parametrized by two angles on the sphere: the azimuthal angle 6{x,t), and the polar angle 
(f){x,t), where 



Iml cos(hsm6, 



Iml sin<?i)sin0, 



\m\ cos 9, 



(7.9) 



and where e [0, 27i), and 9 G [0, vr]. 
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Figure 7.1: The initial data for the magnetization equation (7.3). This initiahzation is obtained 
by allowing the orientation angles of the magnetization vector to vary sinusoidally in space, as in 
Eq. (7.10). Here the wavenumber of the variation is equal to the fundamental wavenumber 2tt/L. 



We carry out numerical simulations of Eqs. (7.3) and (7.5) on a periodic domain [0, L], and outline 



the findings in what follows. Motivated by the change of coordinates (7.9), we choose the initial 
data 

0o(a;) =7r(l + sin(2r7rx/L)), (a;) = ^vr (1 + sin (27rsx/L)) , (7.10) 



where r and s are integers. These data are shown in Fig. 7.1 



Case 1: Numerical simulations of Eq. (7.5). Equation (7.5) is usually solved by explicit or implicit 
finite differences |107] . We solve the equation by these methods, and by the explicit spectral 
method [S]. The accuracy and computational cost is roughly the same in each case, and for 
simplicity, we therefore employ explicit finite differences; it is this method we use throughout 
the chapter. Given the initial conditions (7.10), each component of the magnetization m = 
tends to a constant, the energy 



E 



dx 



dm 



dx 



decays with time, and |m| retains its initial value |mp = 1. After some transience, the decay of 



the energy functional becomes exponential in time. These results are shown in Fig. 7.2 



Case 2: Numerical simulations of Eq. (7.3) with a < (3. Given the smooth initial data (|7.10|), in 
time each component of the magnetization ra = (m^ 



niy, rriz) decays to zero, while the energy 



E = \ I dxm ■ (l - oi^dl) ^ m 
Jq 



tends to a constant value. Given our choice of initial conditions, the energy in fact increases to 
attain this constant value. Again the quantity \m\ stays constant. These results are shown in 
Fig. |7.3[ We find similar results when we set a = 0. 



Case 3: Numerical simulations of Eq. (7.3) with a > (3. Given the smooth initial data (7.10), 
in time each component of the magnetization m = {mx,rny,mz) develops finer and finer scales. 
The development of small scales is driven by the decreasing nature of the energy functional, which 
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Figure 7.2: Numerical simulations of Case (1), the Landau-Lifsliitz-Gilbert equation in the over- 
damped limit. In this case, the magnetization decays to a constant state. Subfigures (a) and (b) 
show the magnetization at times t = 0.03 and t = 0.15 respectively; (c) is the energy functional, 
which exhibits exponential decay after some transience. The final orientation is {4>,6) = {7t,7t/2). 
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Figure 7.3: Numerical simulations of Case (2), the nonlocal Gilbert equation with with a < j3. 
In this case, the energy increases to a constant value, and the magnetization becomes constant. 
Subfigures (a) and (b) show the magnetization at times t = 8 and t = 40; (c) is the energy 
functional. The final orientation is {4>,6) = (vr, 7r/2). 



decreases as power law at late times, and is reflected in snapshots of the power spectrum of 



the magnetization vector, shown in Fig. 7.4 As the system evolves, there is a transfer of large 



Case 


Lengthscales 


Energy 


Outcome as t — > cx) 


Linear Stability 


(1) 


p = o, 

6E/6m = —dim 


Decreasing 


Constant state 


Stable 


(2) 


a < /3 


Increasing 


Constant state 


Stable 


(3) 


a> (3 


Decreasing 


Development of finer 
and finer scales 


Unstable 



Table 7.1: Summary of the forms of Eq. (7.3) studied 



amplitudes to higher wavenumbers. This transfer slows down at late times, suggesting that the rate 
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at which the solution roughens tends to zero, as t ^ oo. The evolution preserves the symmetry of 



the mag netization vector m {x, t) under parity transformations. This is seen by comparing Figs. 7.1 



m 2 

and 7.4 The energy is a decaying function of time, while the quantity \m\ stays constant. We 



find similar results for the case when (3 = 0. 
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Figure 7.4: Numerical simulations of Case (3), the nonlocal Gilbert equation with a > (3. In this 
case, the energy decreases indefinitely, and the magnetization vector develops finer and finer scales. 
Subfigures (a), (b), and (c) show the magnetization at time t = 10000; (d) is the energy functional, 
which decreases in time as a power law at late times. Subfigure (e) shows the power spectrum of 
rrix] the integer index n labels the spatial scales: if /c„ is a wavenumber, then the corresponding 
integer label is n = knL/2'K. 



These results can be explained qualitatively as follows. In Case (1), the energy functional exacts a 
penalty for the formation of gradients. The energy decreases with time and the system evolves into 
a state in which no magnetization gradients are present, that is, a constant state. On the other 
hand, we have demonstrated that in Case (2), when a < (3, the energy increases to a constant 
value. Since in the nonlocal model, the energy functional represents the cost of forming smooth 
spatial structures, an increase in energy produces a smoother magnetization field, a process that 
continues until the magnetization reaches a constant value. Finally, in Case (3), when a > (3, the 
energy functional decreases, and this decrease corresponds to a roughening of the magnetization 



field, as seen in Fig. 7.4, In Sec. 7.4 we shall show that Case (2) is stable to small perturbations 



around a constant state, while Case (3) is unstable. These results are summarized in Table 7.1 



I I 2 Ah. /IK 

conservation of \m\ in Eqs. (7.3) and (7.5) provides a pointwise bound on the magnitude of the 



The solutions of Eqs. (7.3) and (7.5) do not become singular. This is not surprising: the manifest 
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solution, preventing blow-up. Any addition to Eq. (7.3) that breaks this conservation law gives 



rise to the possibility of singular solutions, and it is to this possility that we now turn. 



7.4 Coupled density-magnetization equations 



In this section we study a coupled density-magnetization equation pair that admit singular solu- 
tions. We investigate the linear stability of the equations and examine the conditions for instability. 
We find that the stability or otherwise of a constant state is controlled by that state's magnetiza- 
tion and density, and by the relative magnitude of the problem lengthscales. Using numerical and 
analytical techniques, we investigate the emergence and self-interaction of singular solutions. 



The equations we study are as follows. 



~dt 



d_ 
dx 



dx 6ip ^™ dx 6m 



dm 



d_ 
dx 



d 6E 



where we set 



and, as before. 



5E 



fim = (l-Pidl) ^m 



d 6E 
dx 5m 

-(1- 

5E_ 
5m 



mX \ fJ,mX 



5E 
5m 



(7.11a) 



(7.11b) 



a 



[^-al,dl) ^m 



The order parameter il){x,t) has the interpretation of a particle density, while the orientation 
vector m{x,t) has the interpretation of a local magnetization field. These equations have been 
introduced by Holm, Putkaradze and Tronci in [62] , using a kinetic-theory description. The density 
and the magnetization vector are driven by the velocity 



d 5E 



The velocity advects the ratio \m\/tlj by 



di 



V- 



d 



\m\ 



dx J ip 



d 5E 
dx 5m' 



0. 



(7.12) 



We have the system energy 



E = i / dxm-{l-aldl)-'m-l [ rfx^ (l - aj^^)" V, (7.13) 
Jn Jn 

and, given a nonnegative density, the second term is always nonpositive. This represents an 
energy of attraction, and we therefore expect singularities in the magnetization vector to arise 
from a collapse of the particle density due to the ever-decreasing energy of attraction. There are 
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three lengthscales in the problem that control the time evolution: the ranges am and of the 



potentials in Eq. (7.13), and the smoothening length Pm- 



Linear stability analysis 

We study the linear stability of the constant state (m, ip) = (mo, V'o)- We evaluate the smoothened 
values of this constant solution as follows, 

{l-aldiy^tpo = f{x), 

/ (x) = ^/'o + ^ sinh (a;/a^) + 5 cosh (x/a^) . 
For periodic or infinite boundary conditions, the constants A and B are in fact zero and thus 

(l-a^92)-Vo = V'o, (7.14) 



and similarly hq = {6E/6m] 
(mo,'?/'o) is indeed a solution of Eq. (|7.11[). 



mo 



mo. The result (7.14) guarantees that the constant state 



We study a solution {m,tlj) = (mo + 6m,ijjo + 6ijj), which represents a perturbation away from 
the constant state. By assuming that 6m and 5ijj are initially small in magnitude, we obtain the 
following linearized equations for the perturbation density and magnetization, 

d d"^ -1 d'^ -1 

= -V^0 7T-T7 (1 5ilj + ilJo—{l-al^d'^^ mo -dm, 



dt 



d_ 
dt 



5m = mo 



+ mo X <! mo x 



"(1 - aldlY' 5m - (1 - Pldiy' 5m] }. 



For mo 7^ we may choose two unit vectors ni and n2 such that mo/|mo|, ni and 0,2 form an 
orthonormal triad (that is, we have effected a change of basis). We then study the quantities 5?/^, 
5x, 5^1 and 5^2, where 

5x = n^o ■ S'm, ^^i = '^1 " ^n^^ ^^2 = ^2 ■ 5m. 
We obtain the linear equations 

= -^0^ (1 - '^^ + (1 - "m-^^)"' 5X, 



dt 
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By focussing on a single-mode disturbance with wavenumber k 
equations 



d 
dt 



/ tppk 



|rr^o|^fc2 



V 


















2 = 1,2. 

we obtain the following system of 











1_ 



with eigenvalues 



1 1 / 



/ 



ao = 0, 



l + a2A;2 " 1 + a^F^ 



CTl = 



0-2 = 



l + a^A:2- 



(7.16) 



The eigenvalues are the growth rate of the disturbance {5ip^ 5%, ^^2) |108j . There are two routes 
to instability, when <Ji > 0, or when a2 > 0. The first route leads to an instability when 



1 + 



while the second route leads to instability when 



0-2 > 0, Pm- 



We have plotted the growth rates for the case when ipQ = |moP = 1, and compared the theory with 
numerical simulations. There is excellent agreement at low wavenumbers, although the numerical 
simulations become less accurate at high wavenumbers. This can be remedied by increasing the 



resolution of the simulations. These plots are shown in Figs. 7.5 and 7.6 The growth rates ai^2 
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Figure 7.5: The first route to instability. Subfigure (a) shows the growth rate cxi for am < Pm < c^tp, 
with negativity indicating a stable equilibrium; (b) gives the growth rate cxi for < < Pm-, 
with positivity indicating an unstable equihbrium. We have set |mo| = ^0 = 1- 
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are parabolic in k at small k; ai saturates at large k, while a2 attains a maximum and decays at 
large k. The growth rates can be positive or negative, depending on the initial configuration, and 
on the relationship between the problem lengthscales. In contrast to some standard instabilities of 



pattern formation (e.g. Cahn-Hilliard (Sec. 2.2) or Swift-Hohenberg [109J ), the ai-unstable state 
becomes more unstable at higher wavenumbers (smaller scales), thus preventing the 'freezing-out' 
of the instability by a reduction of the box size |110j . The growth at small scales is limited, 
however, by the saturation in a as ^ oo. Heuristically, this can be explained as follows: at 
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Figure 7.6: The second route to instability. Subfigure (a) shows the growth rate ct2 for am < Pm, 
with negativity indicating a stable equilibrium; (b) gives the growth rate (T2 for am > Pm, with 
positivity indicating an unstable equilibrium. We have set |mo| = V'o = 1- 



higher wavenumber, the disturbance {Sip, Sx,SC,i, 6^,2) gives rise to more and more peaks per unit 
length. This makes merging events increasingly likely, so that peaks combine to form larger peaks, 
enhancing the growth of the disturbance. 



Recall in Sec. 7.3 that the different behaviours of the magnetization equation (7.3) are the result 
of a competition between the lengthscales am and Pm- For am < Pm the initial (large-amplitude) 
disturbance tends to a constant, while for am > Pm the initial disturbance develops finer and 
finer scales. In this section, we have shown that the coupled density-magnetization equations are 
linearly stable when < Pm, while the reverse case is unstable. In contrast to the first route 
to instability, the growth rate a2, if positive, admits a maximum. This is obtained by setting 
CT2 (k) = 0. Then the maximum growth rate occurs at a scale 



A. 



27rk: 



2'K^/amPm- 



Thus, the scale at which the disturbance is most unstable is determined by the geometric mean 
of am and Pm- Given a disturbance {6ip,6x,S^i,6^2) with a range of modes initially present, the 
instability selects the disturbance on the scale Amax- This disturbance develops a large amplitude 
and a singular solution subsequently emerges. It is to this aspect of the problem that we now turn. 
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Singular solutions 

In this section we show that a finite weighted sum of delta functions satisfies the partial differential 



equations (7.11 ). Each delta function has the interpretation of a particle or dumpon, whose weights 



and positions satisfy a finite set of ordinary differential equations. We investigate the two-clumpon 
case analytically and show that the clumpons tend to a state in which they merge, diverge, or are 
separated by a fixed distance. In each case, we determine the final state of the two clumpons' 
relative orientation. 



To verify that singular solutions are possible, let us substitute the ansatz 



M 



M 



tjj {x, t) = ^ai (t) 5{x-Xi (t)) , m{x,t) = ^ hi (t) 5 {x - Xi (t)) 



(7.17) 



i=l 



into the weak form of equations (7.11 ). Here we sum over the different components of the singular 



solution (which we call clumpons). In this section we work on the infinite domain x G (— oo, oo) 



The weak form of the equations is obtained by testing Eqs. (7.11 ) with once-differentiable functions 
f{x) and (/(x). 



dxip (x, t) f (x) 



,,,,,, d 6E d 6E 

dxf {X,t) /i^T^y- + 

' ox dip ox dm 



(7.18a) 



d 
dt 



dxm (x, t) ■ g (x) 



, / . N , M 9 6E d 6E 



+ 



dxg (x) 



m X Ui„ X 



6E 
6m 



(7.18b) 



Substitution of the ansatz (7.17) into the weak equations (7.18) yields the relations 

dxi 



dtti 
~dt 



0, 



dt 



-Vix,) 



dhi 
~dt 



biX ( 



6E 
5m 



ix^) 



ze{l,...,M}, (7.19) 



where V and (^t.^ x {SE/5m)) are obtained from the ansatz (7.17) and are evaluated at Xj. Note 
that the density weights and the magnitude of the weights bi remain constant in time. 

We develop further understanding of the clumpon dynamics by studying the two-clumpon version 
of Eqs. (7.19). Since the weights ai, 02, \bi\, and I62I are constant, two variables suffice to describe 
the interaction: the relative separation X — Xi — X2 of the clumpons, and the cosine of the angle 
between the clumpon magnetizations, cosip = bi ■ b2/|^i||^2|- Using the properties of the kernel 
H (0) = 1, H' (0) = 0, we derive the equations 



dx 
~dt 



MH' (x) - B,H'^^ (x) Hp^ (x) - B,H' (x) y, 2/ = cos ^ 



dy 
dt 



(7.20a) 
(7.20b) 
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where M = 01 + 02, Bi = + |b2p, and B2 = 2|bi||b2| are constants. In this framework, Ha (x) 
is the Green's function of the Helmholtz operator (1 — a^A) , Ha (x) oc e~°'^'. Equations (7.20) 
form a dynamical system whose properties we now investigate using phase-plane analysis [lllj. We 
note first of all that the ||/| > 1 region of the phase plane is forbidden, since the ^/-component of the 




Figure 7.7: The nullcline dx/dt = of the two-clumpon dynamical system with am < cxp- The 
region contained inside the dotted lines ?/ = ±1 gives the allowed values of the dynamical variables 
{x,y). Subfigure (a) shows the case when < C(m- The stable equilibria of the system are 
{x, y) = {±d, 1) and the line x = 0. All initial conditions flow into one of these equilibrium states; 
subfigure (b) shows the case when am < Pm- Initial conditions confined to the line y = 1 flow into 
the fixed point {±d, 1), while all other initial conditions flow into the line x = 0. 

vector field {dx/dt, dy/dt) vanishes at \y\ = 1. The vertical lines x = and x = ±00 are equilibria, 
although their stability will depend on the value of the parameters {am, o^^, jSm, Bi, B2, M) . The 
curve across which dx/dt changes sign is called the nullcline. This is given by 

_ MH'^^ {x) - B,H'^^ {x) Hp^ {x) 

B2H'a^ {X) 

which on the domain x G (—00, 00) takes the form 

a^p O^m 

Several qualitatively different behaviours are possible, depending on the magnitude of the values 
taken by the parameters {am, oixi), Pm, Bi, B2, M). Here we outline four of these behaviour types. 



" 'B~2 



Case 1: The lengthscales are in the relation am < ct^, and (3m < ctm- The nullcline is shown 
in Fig. 7.7 There is flow into the fixed points {x, y) = (±c/, 1), and into the line x = 0, while 
y is a nondecreasing function of time, which follows from Eq. (7.20b). The ultimate state of 
the system is thus x = ±d, cp = (alignment), or x = (merging). In the latter case the 
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Figure 7.8: The nuUcline dx/dt = of the two-clumpon dynamical system with ap < am- The 
region contained inside the dotted hnes y = ±1 gives the allowed values of the dynamical variables 
{x, y). Subfigure (a) shows the case when /3m < am- The lines x = and x = ±oo form the stable 
equilibria of the system. All initial conditions flow into one of these states; subfigure (b) shows the 
case when am < (3m- Initial conditions confined to the line y = 1 flow into the fixed points (0, 1) 
and (±cxD, 1), while all other initial conditions flow into the line x = 0. 



flnal orientation is given by the integral of Eq. (7.20b), 

r /*oo 

tan (I) = tan (^) exp -B^ / dt [Hp^ (x (t)) - H^^ {x (t))] 



<^o = y,(t = 0). (7.21) 



Case 2: The lengthscales are in the relation am < ctip, Cim < f3m, and the nullcline is the 
same as Case (2). All flow along the line y = 1 is directed towards the flxed points {±d, 1). 
All other initial conditions flow into the line x = 0, since y is now a nonincreasing function 
of time. The ultimate state of the system is thus x = ±d, ip = (alignment), or x = 



(merging). In the latter case the flnal orientation is given by the formula (7.21) 



Case 3: The lengthscales are in the relation a^ < a. 
in Fig. [71 



and Pm < ctm- The nullcline is shown 
Inside the region bounded by the line y = and the nullcline, the flow is into 



Case 


am vs. a^ 


am vs. Pm 


Equilibria 


Flow 


(1) 




/3m ^ O^m 


{x,y) = {±d,l); 
X = 0; X = ±oo 


Flow into X = and 

{x,y) = {±d,l) 


(2) 


*^ ^'tp 


O^m ^ (3m 


{x,y) = {±d,l) 
X = 0; X = ±oo 


Flow into X = and 
{x,y) = {±d,l) 


(3) 




(3m ^ Olm 


{x,y) = {±d, 1); 
X = 0; X = ±oo 


Flow into X = 

and X = ±oo 


(4) 






(x,y) = {±d,l); 
X = 0; X = ±oo 


Flow into X = 

and X = ±oo 



Table 7.2: Summary of the distinct phase portraits of Eq. (7.20) studied. 
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the line x = (merging), and the fixed points (±(i, 1) are unstable. The fiow below the line 
y = is towards the line x = 0. Outside of these regions, however, the fiow is into the lines 
X = ±oo, which shows that for a suitable choice of parameters and initial conditions, the 
clumpons can be made to diverge. 

• Case 4-' The lengthscales are in the relation < am and am < Pm, and the nullcline is the 
same as Case (3). The quantity y is a nonincreasing function of time. All fiow along the line 
y = 1 is directed away from the fixed points (±(i, 1) and is into the fixed points (0, 1), or 
(±oo, 1). All other initial conditions fiow into x = 0, although initial conditions that start 
above the curve formed by the nullcline fiow in an arc and eventually reach a fixed point 
(x = 0,2/<0). 



We summarize the cases we have discussed in Table 7.2 Using numerical simulations of Eqs. ( 7.20 ) 



we have verified that Cases (l)-(4) do indeed occur. The list of cases we have considered is not 
exhaustive: depending on the parameters Bi, B2, and M, other phase portraits may arise. Indeed, 
it is clear from Fig. 7.7 that through saddle-node bifurcations, the fixed points {x,y) = {±d,l) 
may disappear, or additional fixed points {x,y) = {±d',—l) may appear. Our analysis shows, 
however, that it is possible to choose a set of parameters {a^, am, Pm, Bi, B2, M) such that two 
clumpons either merge, diverge, or are separated by a fixed distance. 



Numerical Simulations 



To examine the emergence and subsequent interaction of the clumpons, we carry out numerical 
simulations of Eq. (7.11) for a variety of initial conditions. We use an explicit finite-difference 



algorithm with a small amount of artifical diffusion. We solve the following weak form of Eq. (7.11 ) 
obtained by testing Eq. (7.18) with Hp^, 



dip 
~dt 



D. 



art if 



dx"^ 



+ / dyH'{x-y)^{y,t)V{y,t) 



^ = ^artif 1^ + ^ dyH'f,^ (x - y) m, (y, t) V (y, t) 

+ / dyHp^ (x -y)e,- 
Jn 

where ip = Hf^^ * ip and Cj is the unit vector in the i^^ direction. We work on a periodic domain 
Q = [—L/2, L/2], at a resolution of 250 gridpoints; going to higher resolution does not noticeably 
increase the accuracy of the results. 

The first set of initial conditions we study is the following. 



mx Ui^ X — 
V dm J 



m (x, 0) = {sin {AkoX + (px) , sin {Akox + (j)y) , sin {AkoX + (pz)) , 
ip{x,0) = 0.5 + 0.35 cos (2fcox), 



(7.22) 
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where (px, 4>y, and (pz are random phases in the interval [0, 27r], and ko = 2tt/L is the fundamental 
wavenumber. The initial conditions for the magnetization vector are chosen to represent the lack 
of a preferred direction in the problem. The time evolution of equations (7.11) for this set of 



initial conditions is shown in Fig. 7.9 After a short time, the initial data become singular, and 



subsequently, the solution {ijj, m) can be represented as a sum of clumpons. 



M 



M 



ilj{x,t) = ^ai5{x - Xi{t)) , m{x,t) = ^bi{t)6{x - x^{t)) , M = 2. 



1=1 



i=l 



Here M = 2 is the number of clumpons present at the singularity time. This number corresponds 
to the number of maxima in the initial density profile. The forces exerted by each clumpon on the 
other balance because of the effect of the periodic boundary conditions. Indeed, any number of 
equally-spaced, identical, interacting particles arranged on a ring are in equilibrium, although this 
equilibrium is unstable for an attractive force. Thus, at late times, the clumpons are stationary, 
while the magnetization vector fim shows alignment of clumpon magnetizations. 
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Figure 7.9: Evolution of sinusoidally-varying initial conditions for the density and magnetization, 
as in Eq. (7.22). Subfigure (a) shows the evolution of * ip for t E [0,0.15], by which time the 



initial data have formed two clumpons; (b) shows the evolution of fi^- The profiles of jj,y and /i^ are 
similar. Note that the peaks in the density profile correspond to the troughs in the magnetization 
profile. This agrees with the linear stability analysis, wherein disturbances in the density give rise 
to disturbances in the magnetization. 



We gain further understanding of the formation of singular solutions by studying the system 
velocity V just before the onset of the singularity. This is done in Fig. 7.10 Figure 7.10 (a) 



shows the development of the two clumpons from the initial data. Across each density maximum, 
the velocity has the profile V ^ X{t) x, where X{t) > is an increasing function of time. This 
calls to mind the advection problem for the scalar 6 {x,t), studied by Batchelor in the context of 
passive-scalar mixing [22] 

de de 

XqX—, Aq > 0. 



dt 



dx ' 
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Figure 7.10: Evolution of sinusoidally- varying initial conditions for the density and magnetization, 
as in Eq. (7.22). Subfigure (a) shows the system velocity V given in Eq. (7.12), just before the 
singularity time; (b) shows the magnetization at the same time. The density maxima emerge 
at the locations where the convergence of —V (flow into x = and x = ±L/2) occurs, and the 
magnetization develops extrema there. 



Given initial data 9 (x, 0) = 9oe ^^/^o^ the solution evolves in time as 

e(x,t) = ^oe-^M^»^"""*), 
so that gradients are amplified exponentially in time, 



09 
dx 



in a similar manner to the problem studied. 



The evolution of the set of initial conditions (7.22) has therefore demonstrated the following: the 



local velocity V is such that before the onset of the singularity, matter is compressed into regions 
where ip {x, 0) is large, to such an extent that the matter eventually accumulates at isolated points, 
and the singular solution emerges. Moreover, the density maxima, rather than the magnetization 
extrema, drive the formation of singularities. This is not surprising, given that the attractive part 
of the system's energy comes from density variations. 

To highlight the interaction between clumpons, we examine the following set of initial conditions. 



m {x, 0) = tuq = const., ip (x, 0) = 0.5 + 0.35 cos (Skox) 



(7.23) 



where = 27t/L is the fundamental wavenumber. Since this set of initial conditions contains a 
large number of density maxima, we expect a large number of closely-spaced clumpons to emerge. 



and this will illuminate the clumpon interactions. The time evolution of equations (7.11) for this 



set of initial conditions is shown in Fig. 7.11 As before, the solution becomes singular after a 
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Figure 7.11: Evolution of a flat magnetization field and a sinusoidally- varying density, as in 
Eq. (7.23). Subfigure (a) shows the evolution of H^p * ip for t G [0.5,1]; (b) shows the evolu- 
tion of fix- The profiles of Hy and Hz are similar. At t = 0.5, the initial data have formed eight 
equally spaced, identical clumpons, corresponding to the eight density maxima in the initial con- 
figuration. By impulsively shifting the clumpon at x = by a small amount, the equilibrium is 
disrupted and the clumpons merge repeatedly until only one clumpon remains. 



short time, and is subsequently represented by a sum of clumpons, 



M 



M 



ip {x, t) = ^ ai6 {x — Xi (t)) , m (x, t) = ^ 6^ (t) 5 {x — Xi (t)) , M = 8. 
1=1 1=1 

Here M = 8 is the number of clumpons at the singularity time. This number corresponds to the 
number of maxima in the initial density profile. As before, this configuration of equally spaced, 
identical clumpons is an equilibrium state, due to periodic boundary conditions. Therefore, once 
the particle-like state has formed, we impulsively shift the clumpon at x = by a small amount, 
and precipitate the merging of clumpons. The eight clumpons then merge repeatedly until only a 
single clumpon remains. The tendency for the clumpons to merge is explained by the velocity V , 
which changes sign across a clumpon. Thus, if a clumpon is within the range of the force exerted 
by its neighbours, the local velocity, if unbalanced, will advect a given clumpon in the direction of 



one of its neighbours, and the clumpons merge. This process is shown in Fig. 7.13 



7.5 Summary 

We have investigated the nonlocal Gilbert (NG) equation introduced by Holm, Putkaradze, and 
Tronci in [62] using a combination of numerical simulations and simple analytical arguments. The 
NG equation contains two competing lengthscales of nonlocality: there is a lengthscale a associated 
with the range of the interaction potential, and a lengthscale (3 that governs the smoothened 
magnetization vector that appears in the equation. When a < (3 aX\ initial configurations of 
the magnetization tend to a constant value, while for (3 < a the initial configuration of the 
magnetization field develops finer and finer scales. These two effects are in balance when a = /3, 
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Figure 7.12: Evolution of a flat magnetization field and a sinusoidally- varying density, as in 



Eq. (7.23). Subfigure (a) shows the system velocity V given in Eq. (7.12) just before the sin- 



gularity time; (b) gives the magnetization fim, at the same time. The density maxima emerge at 
the locations where the convergence of —V occurs, and the magnetization develops extrema there. 



and the system does not evolve away from its initial state. Furthermore, the NG equation conserves 
the norm of the magnetization vector m, thus providing a pointwise bound on the solution and 
preventing the formation of singular solutions. 

To study the formation of singular solutions, we couple the NG equation to a scalar density 
equation. Associated with the scalar density is a negative energy of attraction that drives the 
formation of singular solutions and breaks the pointwise bound on the m. Three lengthscales of 
nonlocality now enter into the problem: the range of the force associated with the scalar density, 
the range of the force due to the magnetization, and the smoothening length. As before, the 
competition of lengthscales is crucial to the evolution of the system; this is seen in the linear 
stability analysis of the coupled equations, in which the relative magnitude of the lengthscales 
determines the stability or otherwise of a constant state. 

Using numerical simulations, we have demonstrated the emergence of singular solutions from 
smooth initial data, and have explained this behaviour by the negative energy of attraction pro- 
duced by the scalar density. The singular solution consists of a weighted sum of delta functions. 



given in Eq. (7.17), which we interpret as interacting particles or clumpons. The clumpons evolve 
under simple finite-dimensional dynamics. We have shown that a system of two clumpons is gov- 
erned by a two-dimensional dynamical system that has a multiplicity of steady states. Depending 
on the lengthscales of nonlocality and the clumpon weights, the two clumpons can merge, diverge, 
or align and remain separated by a fixed distance. 
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Figure 7.13: Evolution of a flat magnetization field and a sinusoidally- varying density, as in 
Eq. (7.23). Shown is the velocity profile for t G [0.5, 1]; the system velocity is given by Eq. (7.12). 
The velocity —V flows into each density maximum, concentrating matter at isolated points and 
precipitating the formation of eight equally-spaced identical clumpons. On a periodic domain, 
such an arrangement is an equilibrium state, although it is unstable. Thus, by impulsively shifting 
the clumpon at x = by a small amount, we force the clumpons to collapse into larger clumpons, 
until only a single clumpon remains. 



Chapter 8 



Conclusions 



The advective Cahn-Hilliard equation, describing the twin effects of advection and phase separa- 
tion, has been the subject of this report. We have considered cases with and without feedback 
from concentration gradients into the advecting flow field. 

In Ch. [2] we introduced the theory necessary for the report, and one new result: an existence 
theory for the passively-advected Cahn-Hilliard equation. Although an important exercise, such 
existence theories generally tell us no more than the regularity class to which the solution be- 
longs. Therefore, we resorted to numerical simulations in Ch. |3} considering phase separation in 
the presence of a passive, chaotic flow. We characterized the model random-phase sine flow by its 
Lyapunov exponent and correlation length. We quantified the phase separation using the morphol- 
ogy, lengthscales, and probability distribution function of the concentration field. We showed how 
a chaotic flow will not only arrest the growth of the phase-separated domain structure, but will 
lead to remixing of the binary fluid. We highlighted the importance of this result in applications 
where the coarsening tendency of the binary fluid is undesirable. 

Continuing in this vein, in Ch. |4]we defined what it means for a symmetric binary fluid to be well 
mixed and, given a smooth stirring velocity, and sources and sinks, we obtained very precise bounds 
on the level of mixedness achievable by the flow. The question then arises, given a certain source 
term, is there a best flow for mixing a binary fluid? Unlike in the theory of miscible liquids, we 
could not give a complete answer, but noted that the random-phase sine flow of Ch. |3]is efficacious, 
but not necessarily optimal, at mixing. 

In Ch. [5} we switched focus to the active Cahn-Hilliard equation, which is coupled to the Navier- 
Stokes equations. By specializing to the case where the binary fluid forms a thin layer on a 
substrate, it was possible to obtain new, coupled equations, describing the coevolution of the film's 
free-surface height and the fluid concentration, the so-called Stokes Cahn-Hilliard equations. Since 
we were interested in the late-time behaviour of the system, we focussed on the cases where film 
rupture was impossible, by introducing a repulsive Van der Waals force. In that case, we obtained 
existence, regularity, and uniqueness results for solutions to the Stokes Cahn-Hilliard equations 
in two dimensions, and then passed over to numerical simulations in two and three dimensions in 
Ch. [6j The most salient result of this chapter was obtained by applying a spatially- varying surface 



147 



148 



tension across the film. Then, for moderate backreaction strengths, the binary fluid domains 
aligned with the direction of this variation, while for zero backreaction, the domains aligned in the 
perpendicular direction. This surprising result highlights the importance of coupled models, and 
gives a means of controlling phase separation. 

In Ch. [7] we formulated a general theory of aggregation equations that includes as special cases 
the Cahn-Hilliard, Keller-Segel, and (scalar) Holm-Putkaradze equations. Motivated by this 
link between the theory we have developed and the nanoscale physics described by the Holm- 
Putkaradze equation, we studied a coupled vector-scalar system of aggregation equations that 
describes the formation of magnetic nanoparticles. Although this system was introduced by Holm, 
Putkaradze, and Tronci in [22] , its properties had yet to be studied. Therefore, using numerical and 
analytical techniques, we highlighted the competition of lengthscales in a purely magnetic system, 
and demonstrated that singular solutions are never possible in such a model. By passing to the 
coupled vector-scalar model, we demonstrated the emergence of singular solutions from smooth 
initial data. We obtained a set of ordinary differential equations satisified by the discrete elements 
of the singular solution (clumpons). By focussing on the two-clumpon case, we showed rigorously 
that the two clumpons can merge, diverge, or align and remain at a fixed distance. We thus gave 
a full characterization of the Holm-Putkaradze-Tronci model of nanomagnetic particles. 

There are many avenues for future work. The advective Cahn-Hilliard equation could be used 
to give a full picture the breakup and coalescence of droplets with high surface tension, an area 
in which there are many experiments but little theory, while more work could be done in finding 
an optimal flow for mixing the binary fluid. Computational topology |112j could be used to give 
an extra dimension to the characterization of concentration morphology in the presence of flow. 
On the other hand, more experiments on and modelling of the thin-film problem would shed 
further light on this aspect of NSCH dynamics. The relevance of Cahn-Hilliard and thin-film type 
equations to the theory of self-assembly could be probed further. It is the author's hope to work 
on at least some of these problems in the future. 
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Appendix A 



The integral /J^^^ du log cos {u 



We compute the integral 



tt/2 

I — I du log cos (u) . 
'o 



Making the substitution t — cos^ (u) , we have, 
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(A-1) 



where B {x, y) and V (x) are the beta and gamma functions respectively. Using the standard 
identities 

r(i) = i, r(i) = v^, r'(|) = -0F(7 + 2iog2), r'(i) = -7, 

where 7 is the Euler-Mascheroni constant, the integral / becomes 

7 = -|7rlog2. (A-2) 



156 



Appendix B 

The regularized Keller— Segel equation 



The authors Keller and Segel introduced their eponymous equation in a celebrated paper that 
describes the aggregation of the slime mould Dictyostelium discoideum Here we outline the 
derivation of a regularized Keller-Segel equation. This particular model has been considered by 
previous authors |113j . and general discussions of the Keller-Segel model are given in [69 l 1114] . 

The Keller-Segel equation is an equation for the density p{x,t) of a biological population that 
evolves in time through self-advection and diffusion, 



where v = x{p) Vc is the motion of individual elements due to gradients in the concentration field 
c{x,t). This velocity equation can be thought of as a Darcy's Law for a potential — c, although 
in a biological context, it represents the tendency of the population to move towards regions rich 
in the chemical whose concentration is c{x,t). In the situation described by Keller and Segel, the 
chemical c is in fact produced by individual elements in the population. That is, an element of the 
population produces a chemical that attracts its neighbours, which are then advected towards the 
producer. This is encoded in the concentration equation 



where the function g (p, c) is a growth law, here taken to be g (p, c) = —be + ap. Thus, the 
concentration degrades linearly (in time) of its own accord, but increases in the presence of a finite 
population p. In many applications, this concentration field relaxes quickly to a value for which 
dtc = and hence 



dp 
di 



+ V • (vp) = V ■ (DVp) 



dc 
di 



DcAc + g (p, c) , 




or 
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We now examine a Keller-Segel model with a velocity field v = xo{^ ~p/o"o) Vc so that the 
elements of the population do not move once a maximum density do is attained; this takes into 
account the finite size of the individual elements of the population. This approach has been used 
in several studies (see, for example, reference |113j ). and is known to prevent the formation of 
delta-function singularities that would otherwise occur in the density field p{x,t). To find out 
what the appropriate free-energy functional should be, let us work backwards from the evolution 
equation, 



dp 
di 



V • (^DWp - xoP (l - Vc 
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Let us take xo ^i-nd D to be constants. Then, 



dt L V ctq 

= V ■ (M (p) VpKs) , 
Where we have identified the chemical potential 
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(B-1) 



p = 6Fks/Sp = log 



Li-(p/-o)J - f ^ = ^' (^/"°) - 1^' 



and the mobility M {p) = Dp {1 — p/a). The free energy is thus 



Fks = (To / rxW 



Xo 
2D 



d^xp{x,t) J d^yKKs{\x-y\)p{y,t) 



with 



Pks = W {p/ao) - ^Kks *P = W' (p/(7o) - ^c. 



Let us take ip = p/ctq and r = Xot, ^ = D/xq, = Dc/o-Qa and B = b/aoa. Then the model is 
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(B-2) 



which in form is similar to the Cahn-Hilliard equation with double-obstacle potential and variable 
mobility 
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